d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

MFC r271651, r271719, r272138, r272457, r272845, r275476, r275518, r275614,

Browse Source 
r275819, r276176, r278154, r278160, r278339, r279127, r279240, r279491,
    r279493, r279856, r283032, r284423, r284426, r284427, r284428

Merge changes to libm from the past 9 months.  This includes improvements
to the Bessel functions and adds the C99 function lgammal.
This commit is contained in:
tijl 2015-06-25 13:01:10 +00:00
parent 91e9d26e11
commit 95fcbd85e8
43 changed files with 1299 additions and 513 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
lib/msun/Makefile
View File

@ -98,6 +98,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
# If long double != double use these; otherwise, we alias the double versions.
COMMON_SRCS+= e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
e_coshl.c e_fmodl.c e_hypotl.c \
e_lgammal.c e_lgammal_r.c \
e_remainderl.c e_sinhl.c e_sqrtl.c \
invtrig.c k_cosl.c k_sinl.c k_tanl.c \
s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
@ -188,7 +189,8 @@ MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \
ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3
MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \
MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \
lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \
lgamma.3 tgamma.3 lgamma.3 tgammaf.3
MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \
log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \

7
lib/msun/Symbol.map
View File

@ -269,6 +269,7 @@ FBSD_1.3 {
erfl;
expl;
expm1l;
lgammal;
log10l;
log1pl;
log2l;
@ -276,7 +277,11 @@ FBSD_1.3 {
sinhl;
tanhl;
/* Implemented as weak aliases for imprecise versions */
lgammal;
powl;
tgammal;
};
/* First added in 11.0-CURRENT */
FBSD_1.4 {
lgammal_r;
};

330
lib/msun/ld128/e_lgammal_r.c Normal file
View File

@ -0,0 +1,330 @@
/* @(#)e_lgamma_r.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_lgamma_r.c for complete comments.
*
* Converted to long double by Steven G. Kargl.
*/
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
static const volatile double vzero = 0;
static const double
zero= 0,
half= 0.5,
one = 1;
static const long double
pi = 3.14159265358979323846264338327950288e+00L;
/*
* Domain y in [0x1p-119, 0.28], range ~[-1.4065e-36, 1.4065e-36]:
* |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-119.1
*/
static const long double
a0 = 7.72156649015328606065120900824024296e-02L,
a1 = 3.22467033424113218236207583323018498e-01L,
a2 = 6.73523010531980951332460538330282217e-02L,
a3 = 2.05808084277845478790009252803463129e-02L,
a4 = 7.38555102867398526627292839296001626e-03L,
a5 = 2.89051033074152328576829509522483468e-03L,
a6 = 1.19275391170326097618357349881842913e-03L,
a7 = 5.09669524743042462515256340206203019e-04L,
a8 = 2.23154758453578096143609255559576017e-04L,
a9 = 9.94575127818397632126978731542755129e-05L,
a10 = 4.49262367375420471287545895027098145e-05L,
a11 = 2.05072127845117995426519671481628849e-05L,
a12 = 9.43948816959096748454087141447939513e-06L,
a13 = 4.37486780697359330303852050718287419e-06L,
a14 = 2.03920783892362558276037363847651809e-06L,
a15 = 9.55191070057967287877923073200324649e-07L,
a16 = 4.48993286185740853170657139487620560e-07L,
a17 = 2.13107543597620911675316728179563522e-07L,
a18 = 9.70745379855304499867546549551023473e-08L,
a19 = 5.61889970390290257926487734695402075e-08L,
a20 = 6.42739653024130071866684358960960951e-09L,
a21 = 3.34491062143649291746195612991870119e-08L,
a22 = -1.57068547394315223934653011440641472e-08L,
a23 = 1.30812825422415841213733487745200632e-08L;
/*
* Domain x in [tc-0.24, tc+0.28], range ~[-6.3201e-37, 6.3201e-37]:
* |(lgamma(x) - tf) - t(x - tc)| < 2**-120.3.
*/
static const long double
tc = 1.46163214496836234126265954232572133e+00L,
tf = -1.21486290535849608095514557177691584e-01L,
tt = 1.57061739945077675484237837992951704e-36L,
t0 = -1.99238329499314692728655623767019240e-36L,
t1 = -6.08453430711711404116887457663281416e-35L,
t2 = 4.83836122723810585213722380854828904e-01L,
t3 = -1.47587722994530702030955093950668275e-01L,
t4 = 6.46249402389127526561003464202671923e-02L,
t5 = -3.27885410884813055008502586863748063e-02L,
t6 = 1.79706751152103942928638276067164935e-02L,
t7 = -1.03142230366363872751602029672767978e-02L,
t8 = 6.10053602051788840313573150785080958e-03L,
t9 = -3.68456960831637325470641021892968954e-03L,
t10 = 2.25976482322181046611440855340968560e-03L,
t11 = -1.40225144590445082933490395950664961e-03L,
t12 = 8.78232634717681264035014878172485575e-04L,
t13 = -5.54194952796682301220684760591403899e-04L,
t14 = 3.51912956837848209220421213975000298e-04L,
t15 = -2.24653443695947456542669289367055542e-04L,
t16 = 1.44070395420840737695611929680511823e-04L,
t17 = -9.27609865550394140067059487518862512e-05L,
t18 = 5.99347334438437081412945428365433073e-05L,
t19 = -3.88458388854572825603964274134801009e-05L,
t20 = 2.52476631610328129217896436186551043e-05L,
t21 = -1.64508584981658692556994212457518536e-05L,
t22 = 1.07434583475987007495523340296173839e-05L,
t23 = -7.03070407519397260929482550448878399e-06L,
t24 = 4.60968590693753579648385629003100469e-06L,
t25 = -3.02765473778832036018438676945512661e-06L,
t26 = 1.99238771545503819972741288511303401e-06L,
t27 = -1.31281299822614084861868817951788579e-06L,
t28 = 8.60844432267399655055574642052370223e-07L,
t29 = -5.64535486432397413273248363550536374e-07L,
t30 = 3.99357783676275660934903139592727737e-07L,
t31 = -2.95849029193433121795495215869311610e-07L,
t32 = 1.37790144435073124976696250804940384e-07L;
/*
* Domain y in [-0.1, 0.232], range ~[-1.4046e-37, 1.4181e-37]:
* |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-122.8
*/
static const long double
u0 = -7.72156649015328606065120900824024311e-02L,
u1 = 4.24082772271938167430983113242482656e-01L,
u2 = 2.96194003481457101058321977413332171e+00L,
u3 = 6.49503267711258043997790983071543710e+00L,
u4 = 7.40090051288150177152835698948644483e+00L,
u5 = 4.94698036296756044610805900340723464e+00L,
u6 = 2.00194224610796294762469550684947768e+00L,
u7 = 4.82073087750608895996915051568834949e-01L,
u8 = 6.46694052280506568192333848437585427e-02L,
u9 = 4.17685526755100259316625348933108810e-03L,
u10 = 9.06361003550314327144119307810053410e-05L,
v1 = 5.15937098592887275994320496999951947e+00L,
v2 = 1.14068418766251486777604403304717558e+01L,
v3 = 1.41164839437524744055723871839748489e+01L,
v4 = 1.07170702656179582805791063277960532e+01L,
v5 = 5.14448694179047879915042998453632434e+00L,
v6 = 1.55210088094585540637493826431170289e+00L,
v7 = 2.82975732849424562719893657416365673e-01L,
v8 = 2.86424622754753198010525786005443539e-02L,
v9 = 1.35364253570403771005922441442688978e-03L,
v10 = 1.91514173702398375346658943749580666e-05L,
v11 = -3.25364686890242327944584691466034268e-08L;
/*
* Domain x in (2, 3], range ~[-1.3341e-36, 1.3536e-36]:
* |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-120.1
* with y = x - 2.
*/
static const long double
s0 = -7.72156649015328606065120900824024297e-02L,
s1 = 1.23221687850916448903914170805852253e-01L,
s2 = 5.43673188699937239808255378293820020e-01L,
s3 = 6.31998137119005233383666791176301800e-01L,
s4 = 3.75885340179479850993811501596213763e-01L,
s5 = 1.31572908743275052623410195011261575e-01L,
s6 = 2.82528453299138685507186287149699749e-02L,
s7 = 3.70262021550340817867688714880797019e-03L,
s8 = 2.83374000312371199625774129290973648e-04L,
s9 = 1.15091830239148290758883505582343691e-05L,
s10 = 2.04203474281493971326506384646692446e-07L,
s11 = 9.79544198078992058548607407635645763e-10L,
r1 = 2.58037466655605285937112832039537492e+00L,
r2 = 2.86289413392776399262513849911531180e+00L,
r3 = 1.78691044735267497452847829579514367e+00L,
r4 = 6.89400381446725342846854215600008055e-01L,
r5 = 1.70135865462567955867134197595365343e-01L,
r6 = 2.68794816183964420375498986152766763e-02L,
r7 = 2.64617234244861832870088893332006679e-03L,
r8 = 1.52881761239180800640068128681725702e-04L,
r9 = 4.63264813762296029824851351257638558e-06L,
r10 = 5.89461519146957343083848967333671142e-08L,
r11 = 1.79027678176582527798327441636552968e-10L;
/*
* Domain z in [8, 0x1p70], range ~[-9.8214e-35, 9.8214e-35]:
* |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-113.0
*/
static const long double
w0 = 4.18938533204672741780329736405617738e-01L,
w1 = 8.33333333333333333333333333332852026e-02L,
w2 = -2.77777777777777777777777727810123528e-03L,
w3 = 7.93650793650793650791708939493907380e-04L,
w4 = -5.95238095238095234390450004444370959e-04L,
w5 = 8.41750841750837633887817658848845695e-04L,
w6 = -1.91752691752396849943172337347259743e-03L,
w7 = 6.41025640880333069429106541459015557e-03L,
w8 = -2.95506530801732133437990433080327074e-02L,
w9 = 1.79644237328444101596766586979576927e-01L,
w10 = -1.39240539108367641920172649259736394e+00L,
w11 = 1.33987701479007233325288857758641761e+01L,
w12 = -1.56363596431084279780966590116006255e+02L,
w13 = 2.14830978044410267201172332952040777e+03L,
w14 = -3.28636067474227378352761516589092334e+04L,
w15 = 5.06201257747865138432663574251462485e+05L,
w16 = -6.79720123352023636706247599728048344e+06L,
w17 = 6.57556601705472106989497289465949255e+07L,
w18 = -3.26229058141181783534257632389415580e+08L;
static long double
sin_pil(long double x)
{
volatile long double vz;
long double y,z;
uint64_t lx, n;
uint16_t hx;
y = -x;
vz = y+0x1.p112;
z = vz-0x1.p112;
if (z == y)
return zero;
vz = y+0x1.p110;
EXTRACT_LDBL128_WORDS(hx,lx,n,vz);
z = vz-0x1.p110;
if (z > y) {
z -= 0.25;
n--;
}
n &= 7;
y = y - z + n * 0.25;
switch (n) {
case 0: y = __kernel_sinl(pi*y,zero,0); break;
case 1:
case 2: y = __kernel_cosl(pi*(0.5-y),zero); break;
case 3:
case 4: y = __kernel_sinl(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break;
default: y = __kernel_sinl(pi*(y-2.0),zero,0); break;
}
return -y;
}
long double
lgammal_r(long double x, int *signgamp)
{
long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
uint64_t llx,lx;
int i;
uint16_t hx,ix;
EXTRACT_LDBL128_WORDS(hx,lx,llx,x);
/* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fff;
if(ix==0x7fff) return x*x;
/* purge +-0 and tiny arguments */
*signgamp = 1-2*(hx>>15);
if(ix<0x3fff-116) { /* |x|<2**-(p+3), return -log(|x|) */
if((ix|lx|llx)==0)
return one/vzero;
return -logl(fabsl(x));
}
/* purge negative integers and start evaluation for other x < 0 */
if(hx&0x8000) {
*signgamp = 1;
if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */
return one/vzero;
t = sin_pil(x);
if(t==zero) return one/vzero;
nadj = logl(pi/fabsl(t*x));
if(t<zero) *signgamp = -1;
x = -x;
}
/* purge 1 and 2 */
if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0;
/* for x < 2.0 */
else if(ix<0x4000) {
if(x<=8.9999961853027344e-01) {
r = -logl(x);
if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
else {y = x; i=2;}
} else {
r = 0;
if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
else {y=x-1;i=2;}
}
switch(i) {
case 0:
z = y*y;
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+
z*(a18+z*(a20+z*a22))))))))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+
z*(a17+z*(a19+z*(a21+z*a23)))))))))));
p = y*p1+p2;
r += p-y/2; break;
case 1:
p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+
y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+
y*(t31+y*t32))))))))))))))))))))))))))))));
r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+
y*(u8+y*(u9+y*u10))))))))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+
y*(v8+y*(v9+y*(v10+y*v11))))))))));
r += p1/p2-y/2;
}
}
/* x < 8.0 */
else if(ix<0x4002) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+
y*(s9+y*(s10+y*s11)))))))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+
y*(r9+y*(r10+y*r11))))))))));
r = y/2+p/q;
z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+6); /* FALLTHRU */
case 6: z *= (y+5); /* FALLTHRU */
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += logl(z); break;
}
/* 8.0 <= x < 2**(p+3) */
} else if (ix<0x3fff+116) {
t = logl(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*(w8+
y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+
y*(w17+y*w18)))))))))))))))));
r = (x-half)*(t-one)+w;
/* 2**(p+3) <= x <= inf */
} else
r = x*(logl(x)-1);
if(hx&0x8000) r = nadj - r;
return r;
}

2
lib/msun/ld128/k_expl.h
View File

@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int expt)
scale2 = 1;
SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
return (cpackl(cos(y) * exp_x * scale1 * scale2,
return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
sinl(y) * exp_x * scale1 * scale2));
}
#endif /* _COMPLEX_H */

358
lib/msun/ld80/e_lgammal_r.c Normal file
View File

@ -0,0 +1,358 @@
/* @(#)e_lgamma_r.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_lgamma_r.c for complete comments.
*
* Converted to long double by Steven G. Kargl.
*/
#ifdef __i386__
#include <ieeefp.h>
#endif
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
static const volatile double vzero = 0;
static const double
zero= 0,
half= 0.5,
one = 1;
static const union IEEEl2bits
piu = LD80C(0xc90fdaa22168c235, 1, 3.14159265358979323851e+00L);
#define pi (piu.e)
/*
* Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]:
* |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9
*/
static const union IEEEl2bits
a0u = LD80C(0x9e233f1bed863d26, -4, 7.72156649015328606028e-02L),
a1u = LD80C(0xa51a6625307d3249, -2, 3.22467033424113218889e-01L),
a2u = LD80C(0x89f000d2abafda8c, -4, 6.73523010531979398946e-02L),
a3u = LD80C(0xa8991563eca75f26, -6, 2.05808084277991211934e-02L),
a4u = LD80C(0xf2027e10634ce6b6, -8, 7.38555102796070454026e-03L),
a5u = LD80C(0xbd6eb76dd22187f4, -9, 2.89051035162703932972e-03L),
a6u = LD80C(0x9c562ab05e0458ed, -10, 1.19275351624639999297e-03L),
a7u = LD80C(0x859baed93ee48e46, -11, 5.09674593842117925320e-04L),
a8u = LD80C(0xe9f28a4432949af2, -13, 2.23109648015769155122e-04L),
a9u = LD80C(0xd12ad0d9b93c6bb0, -14, 9.97387167479808509830e-05L),
a10u= LD80C(0xb7522643c78a219b, -15, 4.37071076331030136818e-05L),
a11u= LD80C(0xca024dcdece2cb79, -16, 2.40813493372040143061e-05L),
a12u= LD80C(0xbb90fb6968ebdbf9, -19, 2.79495621083634031729e-06L),
a13u= LD80C(0xba1c9ffeeae07b37, -17, 1.10931287015513924136e-05L);
#define a0 (a0u.e)
#define a1 (a1u.e)
#define a2 (a2u.e)
#define a3 (a3u.e)
#define a4 (a4u.e)
#define a5 (a5u.e)
#define a6 (a6u.e)
#define a7 (a7u.e)
#define a8 (a8u.e)
#define a9 (a9u.e)
#define a10 (a10u.e)
#define a11 (a11u.e)
#define a12 (a12u.e)
#define a13 (a13u.e)
/*
* Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]:
* |(lgamma(x) - tf) - t(x - tc)| < 2**-70.5
*/
static const union IEEEl2bits
tcu = LD80C(0xbb16c31ab5f1fb71, 0, 1.46163214496836234128e+00L),
tfu = LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L),
ttu = LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L),
t0u = LD80C(0x80b9406556a62a6b, -68, 3.40728634996055147231e-21L),
t1u = LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L),
t2u = LD80C(0xf7b95e4771c55d51, -2, 4.83836122723810583532e-01L),
t3u = LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L),
t4u = LD80C(0x845a14a6a81dc94b, -4, 6.46249402389135358063e-02L),
t5u = LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L),
t6u = LD80C(0x93373cbd00297438, -6, 1.79706751150707171293e-02L),
t7u = LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L),
t8u = LD80C(0xc7e7015ff4bc45af, -8, 6.10053603296546099193e-03L),
t9u = LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L),
t10u = LD80C(0x94188d58f12e5e9f, -9, 2.25976420273774583089e-03L),
t11u = LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L),
t12u = LD80C(0xe63a671e6704ea4d, -11, 8.78250640744776944887e-04L),
t13u = LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L),
t14u = LD80C(0xb858f5bdb79276fe, -12, 3.51614951536825927370e-04L),
t15u = LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L),
t16u = LD80C(0x99aeabb0d67ba835, -13, 1.46562869351659194136e-04L),
t17u = LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L),
t18u = LD80C(0xe24cb1e3b0474775, -15, 5.39540265505221957652e-05L);
#define tc (tcu.e)
#define tf (tfu.e)
#define tt (ttu.e)
#define t0 (t0u.e)
#define t1 (t1u.e)
#define t2 (t2u.e)
#define t3 (t3u.e)
#define t4 (t4u.e)
#define t5 (t5u.e)
#define t6 (t6u.e)
#define t7 (t7u.e)
#define t8 (t8u.e)
#define t9 (t9u.e)
#define t10 (t10u.e)
#define t11 (t11u.e)
#define t12 (t12u.e)
#define t13 (t13u.e)
#define t14 (t14u.e)
#define t15 (t15u.e)
#define t16 (t16u.e)
#define t17 (t17u.e)
#define t18 (t18u.e)
/*
* Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]:
* |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2
*/
static const union IEEEl2bits
u0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
u1u = LD80C(0x98280ee45e4ddd3d, -1, 5.94361239198682739769e-01L),
u2u = LD80C(0xe330c8ead4130733, 0, 1.77492629495841234275e+00L),
u3u = LD80C(0xd4a213f1a002ec52, 0, 1.66119622514818078064e+00L),
u4u = LD80C(0xa5a9ca6f5bc62163, -1, 6.47122051417476492989e-01L),
u5u = LD80C(0xc980e49cd5b019e6, -4, 9.83903751718671509455e-02L),
u6u = LD80C(0xff636a8bdce7025b, -9, 3.89691687802305743450e-03L),
v1u = LD80C(0xbd109c533a19fbf5, 1, 2.95413883330948556544e+00L),
v2u = LD80C(0xd295cbf96f31f099, 1, 3.29039286955665403176e+00L),
v3u = LD80C(0xdab8bcfee40496cb, 0, 1.70876276441416471410e+00L),
v4u = LD80C(0xd2f2dc3638567e9f, -2, 4.12009126299534668571e-01L),
v5u = LD80C(0xa07d9b0851070f41, -5, 3.91822868305682491442e-02L),
v6u = LD80C(0xe3cd8318f7adb2c4, -11, 8.68998648222144351114e-04L);
#define u0 (u0u.e)
#define u1 (u1u.e)
#define u2 (u2u.e)
#define u3 (u3u.e)
#define u4 (u4u.e)
#define u5 (u5u.e)
#define u6 (u6u.e)
#define v1 (v1u.e)
#define v2 (v2u.e)
#define v3 (v3u.e)
#define v4 (v4u.e)
#define v5 (v5u.e)
#define v6 (v6u.e)
/*
* Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]:
* |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3
* with y = x - 2.
*/
static const union IEEEl2bits
s0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
s1u = LD80C(0xd3ff0dcc7fa91f94, -3, 2.07027640921219389860e-01L),
s2u = LD80C(0xb2bb62782478ef31, -2, 3.49085881391362090549e-01L),
s3u = LD80C(0xb49f7438c4611a74, -3, 1.76389518704213357954e-01L),
s4u = LD80C(0x9a957008fa27ecf9, -5, 3.77401710862930008071e-02L),
s5u = LD80C(0xda9b389a6ca7a7ac, -9, 3.33566791452943399399e-03L),
s6u = LD80C(0xbc7a2263faf59c14, -14, 8.98728786745638844395e-05L),
r1u = LD80C(0xbf5cff5b11477d4d, 0, 1.49502555796294337722e+00L),
r2u = LD80C(0xd9aec89de08e3da6, -1, 8.50323236984473285866e-01L),
r3u = LD80C(0xeab7ae5057c443f9, -3, 2.29216312078225806131e-01L),
r4u = LD80C(0xf29707d9bd2b1e37, -6, 2.96130326586640089145e-02L),
r5u = LD80C(0xd376c2f09736c5a3, -10, 1.61334161411590662495e-03L),
r6u = LD80C(0xc985983d0cd34e3d, -16, 2.40232770710953450636e-05L),
r7u = LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L);
#define s0 (s0u.e)
#define s1 (s1u.e)
#define s2 (s2u.e)
#define s3 (s3u.e)
#define s4 (s4u.e)
#define s5 (s5u.e)
#define s6 (s6u.e)
#define r1 (r1u.e)
#define r2 (r2u.e)
#define r3 (r3u.e)
#define r4 (r4u.e)
#define r5 (r5u.e)
#define r6 (r6u.e)
#define r7 (r7u.e)
/*
* Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]:
* |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7
*/
static const union IEEEl2bits
w0u = LD80C(0xd67f1c864beb4a69, -2, 4.18938533204672741776e-01L),
w1u = LD80C(0xaaaaaaaaaaaaaaa1, -4, 8.33333333333333332678e-02L),
w2u = LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L),
w3u = LD80C(0xd00d00cf58aede4c, -11, 7.93650793490637233668e-04L),
w4u = LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L),
w5u = LD80C(0xdca7cadc5baa517b, -11, 8.41733700408000822962e-04L),
w6u = LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L),
w7u = LD80C(0xcbd5101bb58d1f2b, -8, 6.22046743903262649294e-03L),
w8u = LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L);
#define w0 (w0u.e)
#define w1 (w1u.e)
#define w2 (w2u.e)
#define w3 (w3u.e)
#define w4 (w4u.e)
#define w5 (w5u.e)
#define w6 (w6u.e)
#define w7 (w7u.e)
#define w8 (w8u.e)
static long double
sin_pil(long double x)
{
volatile long double vz;
long double y,z;
uint64_t n;
uint16_t hx;
y = -x;
vz = y+0x1p63;
z = vz-0x1p63;
if (z == y)
return zero;
vz = y+0x1p61;
EXTRACT_LDBL80_WORDS(hx,n,vz);
z = vz-0x1p61;
if (z > y) {
z -= 0.25; /* adjust to round down */
n--;
}
n &= 7; /* octant of y mod 2 */
y = y - z + n * 0.25; /* y mod 2 */
switch (n) {
case 0: y = __kernel_sinl(pi*y,zero,0); break;
case 1:
case 2: y = __kernel_cosl(pi*(0.5-y),zero); break;
case 3:
case 4: y = __kernel_sinl(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break;
default: y = __kernel_sinl(pi*(y-2.0),zero,0); break;
}
return -y;
}
long double
lgammal_r(long double x, int *signgamp)
{
long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
uint64_t lx;
int i;
uint16_t hx,ix;
EXTRACT_LDBL80_WORDS(hx,lx,x);
/* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fff;
if(ix==0x7fff) return x*x;
ENTERI();
/* purge +-0 and tiny arguments */
*signgamp = 1-2*(hx>>15);
if(ix<0x3fff-67) { /* |x|<2**-(p+3), return -log(|x|) */
if((ix|lx)==0)
RETURNI(one/vzero);
RETURNI(-logl(fabsl(x)));
}
/* purge negative integers and start evaluation for other x < 0 */
if(hx&0x8000) {
*signgamp = 1;
if(ix>=0x3fff+63) /* |x|>=2**(p-1), must be -integer */
RETURNI(one/vzero);
t = sin_pil(x);
if(t==zero) RETURNI(one/vzero); /* -integer */
nadj = logl(pi/fabsl(t*x));
if(t<zero) *signgamp = -1;
x = -x;
}
/* purge 1 and 2 */
if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0;
/* for x < 2.0 */
else if(ix<0x4000) {
/*
* XXX Supposedly, one can use the following information to replace the
* XXX FP rational expressions. A similar approach is appropriate
* XXX for ld128, but one (may need?) needs to consider llx, too.
*
* 8.9999961853027344e-01 3ffe e666600000000000
* 7.3159980773925781e-01 3ffe bb4a200000000000
* 2.3163998126983643e-01 3ffc ed33080000000000
* 1.7316312789916992e+00 3fff dda6180000000000
* 1.2316322326660156e+00 3fff 9da6200000000000
*/
if(x<8.9999961853027344e-01) {
r = -logl(x);
if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
else {y = x; i=2;}
} else {
r = 0;
if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
else {y=x-1;i=2;}
}
switch(i) {
case 0:
z = y*y;
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12)))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13))))));
p = y*p1+p2;
r += p-y/2; break;
case 1:
p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
y*(t17+y*t18))))))))))))))));
r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6))))));
p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6)))));
r += p1/p2-y/2;
}
}
/* x < 8.0 */
else if(ix<0x4002) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7))))));
r = y/2+p/q;
z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+6); /* FALLTHRU */
case 6: z *= (y+5); /* FALLTHRU */
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += logl(z); break;
}
/* 8.0 <= x < 2**(p+3) */
} else if (ix<0x3fff+67) {
t = logl(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8)))))));
r = (x-half)*(t-one)+w;
/* 2**(p+3) <= x <= inf */
} else
r = x*(logl(x)-1);
if(hx&0x8000) r = nadj - r;
RETURNI(r);
}

2
lib/msun/ld80/k_expl.h
View File

@ -299,7 +299,7 @@ __ldexp_cexpl(long double complex z, int expt)
scale2 = 1;
SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
return (cpackl(cos(y) * exp_x * scale1 * scale2,
return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
sinl(y) * exp_x * scale1 * scale2));
}
#endif /* _COMPLEX_H */

45
lib/msun/man/j0.3
View File

@ -28,7 +28,7 @@
.\" from: @(#)j0.3 6.7 (Berkeley) 4/19/91
.\" $FreeBSD$
.\"
.Dd February 18, 2008
.Dd March 10, 2015
.Dt J0 3
.Os
.Sh NAME
@ -77,24 +77,17 @@
The functions
.Fn j0 ,
.Fn j0f ,
.Fn j1
.Fn j1 ,
and
.Fn j1f
compute the
.Em Bessel function of the first kind of the order
0 and the
.Em order
1, respectively,
for the
real value
compute the Bessel function of the first kind of orders
0 and 1 for the real value
.Fa x ;
the functions
.Fn jn
and
.Fn jnf
compute the
.Em Bessel function of the first kind of the integer
.Em order
compute the Bessel function of the first kind of the integer order
.Fa n
for the real value
.Fa x .
@ -105,13 +98,8 @@ The functions
.Fn y1 ,
and
.Fn y1f
compute the linearly independent
.Em Bessel function of the second kind of the order
0 and the
.Em order
1, respectively,
for the
positive
compute the linearly independent Bessel function of the second kind
of orders 0 and 1 for the positive
.Em real
value
.Fa x ;
@ -119,9 +107,7 @@ the functions
.Fn yn
and
.Fn ynf
compute the
.Em Bessel function of the second kind for the integer
.Em order
compute the Bessel function of the second kind for the integer order
.Fa n
for the positive
.Em real
@ -141,11 +127,20 @@ and
.Fn ynf .
If
.Fa x
is negative, these routines will generate an invalid exception and
return \*(Na.
is negative, including -\*(If, these routines will generate an invalid
exception and return \*(Na.
If
.Fa x
is 0 or a sufficiently small positive number, these routines
is \*(Pm0, these routines
will generate a divide-by-zero exception and return -\*(If.
If
.Fa x
is a sufficiently small positive number, then
.Fn y1 ,
.Fn y1f ,
.Fn yn ,
and
.Fn ynf
will generate an overflow exception and return -\*(If.
.Sh SEE ALSO
.Xr math 3

37
lib/msun/man/lgamma.3
View File

@ -28,7 +28,7 @@
.\" from: @(#)lgamma.3 6.6 (Berkeley) 12/3/92
.\" $FreeBSD$
.\"
.Dd January 14, 2005
.Dd September 12, 2014
.Dt LGAMMA 3
.Os
.Sh NAME
@ -36,6 +36,8 @@
.Nm lgamma_r ,
.Nm lgammaf ,
.Nm lgammaf_r ,
.Nm lgammal ,
.Nm lgammal_r ,
.Nm gamma ,
.Nm gamma_r ,
.Nm gammaf ,
@ -58,6 +60,10 @@
.Fn lgammaf "float x"
.Ft float
.Fn lgammaf_r "float x" "int *signgamp"
.Ft "long double"
.Fn lgammal "long double x"
.Ft "long double"
.Fn lgammal_r "long double x" "int *signgamp"
.Ft double
.Fn gamma "double x"
.Ft double
@ -66,14 +72,15 @@
.Fn gammaf "float x"
.Ft float
.Fn gammaf_r "float x" "int *signgamp"
.Ft double
.Ft "long double"
.Fn tgamma "double x"
.Ft float
.Fn tgammaf "float x"
.Sh DESCRIPTION
.Fn lgamma x
.Fn lgamma x ,
.Fn lgammaf x ,
and
.Fn lgammaf x
.Fn lgammal x
.if t \{\
return ln\||\(*G(x)| where
.Bd -unfilled -offset indent
@ -87,13 +94,15 @@ The external integer
.Fa signgam
returns the sign of \(*G(x).
.Pp
.Fn lgamma_r x signgamp
.Fn lgamma_r x signgamp ,
.Fn lgammaf_r x signgamp ,
and
.Fn lgammaf_r x signgamp
.Fn lgammal_r x signgamp
provide the same functionality as
.Fn lgamma x
.Fn lgamma x ,
.Fn lgammaf x ,
and
.Fn lgammaf x
.Fn lgammal x ,
but the caller must provide an integer to store the sign of \(*G(x).
.Pp
The
@ -115,6 +124,7 @@ are deprecated aliases for
and
.Fn lgammaf_r ,
respectively.
.Sh IDIOSYNCRASIES
Do not use the expression
.Dq Li signgam\(**exp(lgamma(x))
@ -139,14 +149,18 @@ Exponentiation of
will lose up to 10 significant bits.
.Sh RETURN VALUES
.Fn gamma ,
.Fn gamma_r ,
.Fn gammaf ,
.Fn gammal ,
.Fn gamma_r ,
.Fn gammaf_r ,
.Fn gammal_r ,
.Fn lgamma ,
.Fn lgamma_r ,
.Fn lgammaf ,
.Fn lgammal ,
.Fn lgamma_r ,
.Fn lgammaf_r ,
and
.Fn lgammaf_r
.Fn lgammal_r
return appropriate values unless an argument is out of range.
Overflow will occur for sufficiently large positive values, and
non-positive integers.
@ -159,6 +173,7 @@ will underflow.
The
.Fn lgamma ,
.Fn lgammaf ,
.Fn lgammal ,
.Fn tgamma ,
and
.Fn tgammaf

64
lib/msun/src/catrig.c
View File

@ -286,19 +286,19 @@ casinh(double complex z)
if (isnan(x) || isnan(y)) {
/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
if (isinf(x))
return (cpack(x, y + y));
return (CMPLX(x, y + y));
/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
if (isinf(y))
return (cpack(y, x + x));
return (CMPLX(y, x + x));
/* casinh(NaN + I*0) = NaN + I*0 */
if (y == 0)
return (cpack(x + x, y));
return (CMPLX(x + x, y));
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -307,7 +307,7 @@ casinh(double complex z)
w = clog_for_large_values(z) + m_ln2;
else
w = clog_for_large_values(-z) + m_ln2;
return (cpack(copysign(creal(w), x), copysign(cimag(w), y)));
return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
}
/* Avoid spuriously raising inexact for z = 0. */
@ -325,7 +325,7 @@ casinh(double complex z)
ry = asin(B);
else
ry = atan2(new_y, sqrt_A2my2);
return (cpack(copysign(rx, x), copysign(ry, y)));
return (CMPLX(copysign(rx, x), copysign(ry, y)));
}
/*
@ -335,9 +335,9 @@ casinh(double complex z)
double complex
casin(double complex z)
{
double complex w = casinh(cpack(cimag(z), creal(z)));
double complex w = casinh(CMPLX(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
return (CMPLX(cimag(w), creal(w)));
}
/*
@ -370,19 +370,19 @@ cacos(double complex z)
if (isnan(x) || isnan(y)) {
/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
if (isinf(x))
return (cpack(y + y, -INFINITY));
return (CMPLX(y + y, -INFINITY));
/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
if (isinf(y))
return (cpack(x + x, -y));
return (CMPLX(x + x, -y));
/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
if (x == 0)
return (cpack(pio2_hi + pio2_lo, y + y));
return (CMPLX(pio2_hi + pio2_lo, y + y));
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -392,18 +392,18 @@ cacos(double complex z)
ry = creal(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
return (CMPLX(rx, ry));
}
/* Avoid spuriously raising inexact for z = 1. */
if (x == 1 && y == 0)
return (cpack(0, -y));
return (CMPLX(0, -y));
/* All remaining cases are inexact. */
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (cpack(pio2_hi - (x - pio2_lo), -y));
return (CMPLX(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
@ -419,7 +419,7 @@ cacos(double complex z)
}
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
return (CMPLX(rx, ry));
}
/*
@ -437,15 +437,15 @@ cacosh(double complex z)
ry = cimag(w);
/* cacosh(NaN + I*NaN) = NaN + I*NaN */
if (isnan(rx) && isnan(ry))
return (cpack(ry, rx));
return (CMPLX(ry, rx));
/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
if (isnan(rx))
return (cpack(fabs(ry), rx));
return (CMPLX(fabs(ry), rx));
/* cacosh(0 + I*NaN) = NaN + I*NaN */
if (isnan(ry))
return (cpack(ry, ry));
return (cpack(fabs(ry), copysign(rx, cimag(z))));
return (CMPLX(ry, ry));
return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
}
/*
@ -475,16 +475,16 @@ clog_for_large_values(double complex z)
* this method is still poor since it is uneccessarily slow.
*/
if (ax > DBL_MAX / 2)
return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
/*
* Avoid overflow when x or y is large. Avoid underflow when x or
* y is small.
*/
if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
return (cpack(log(hypot(x, y)), atan2(y, x)));
return (CMPLX(log(hypot(x, y)), atan2(y, x)));
return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
}
/*
@ -575,30 +575,30 @@ catanh(double complex z)
/* This helps handle many cases. */
if (y == 0 && ax <= 1)
return (cpack(atanh(x), y));
return (CMPLX(atanh(x), y));
/* To ensure the same accuracy as atan(), and to filter out z = 0. */
if (x == 0)
return (cpack(x, atan(y)));
return (CMPLX(x, atan(y)));
if (isnan(x) || isnan(y)) {
/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
if (isinf(x))
return (cpack(copysign(0, x), y + y));
return (CMPLX(copysign(0, x), y + y));
/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
if (isinf(y))
return (cpack(copysign(0, x),
return (CMPLX(copysign(0, x),
copysign(pio2_hi + pio2_lo, y)));
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
return (cpack(real_part_reciprocal(x, y),
return (CMPLX(real_part_reciprocal(x, y),
copysign(pio2_hi + pio2_lo, y)));
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@ -623,7 +623,7 @@ catanh(double complex z)
else
ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpack(copysign(rx, x), copysign(ry, y)));
return (CMPLX(copysign(rx, x), copysign(ry, y)));
}
/*
@ -633,7 +633,7 @@ catanh(double complex z)
double complex
catan(double complex z)
{
double complex w = catanh(cpack(cimag(z), creal(z)));
double complex w = catanh(CMPLX(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
return (CMPLX(cimag(w), creal(w)));
}

64
lib/msun/src/catrigf.c
View File

@ -156,12 +156,12 @@ casinhf(float complex z)
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(x, y + y));
return (CMPLXF(x, y + y));
if (isinf(y))
return (cpackf(y, x + x));
return (CMPLXF(y, x + x));
if (y == 0)
return (cpackf(x + x, y));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(x + x, y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -169,7 +169,7 @@ casinhf(float complex z)
w = clog_for_large_values(z) + m_ln2;
else
w = clog_for_large_values(-z) + m_ln2;
return (cpackf(copysignf(crealf(w), x),
return (CMPLXF(copysignf(crealf(w), x),
copysignf(cimagf(w), y)));
}
@ -186,15 +186,15 @@ casinhf(float complex z)
ry = asinf(B);
else
ry = atan2f(new_y, sqrt_A2my2);
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
}
float complex
casinf(float complex z)
{
float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
return (CMPLXF(cimagf(w), crealf(w)));
}
float complex
@ -214,12 +214,12 @@ cacosf(float complex z)
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(y + y, -INFINITY));
return (CMPLXF(y + y, -INFINITY));
if (isinf(y))
return (cpackf(x + x, -y));
return (CMPLXF(x + x, -y));
if (x == 0)
return (cpackf(pio2_hi + pio2_lo, y + y));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(pio2_hi + pio2_lo, y + y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -228,16 +228,16 @@ cacosf(float complex z)
ry = crealf(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpackf(rx, ry));
return (CMPLXF(rx, ry));
}
if (x == 1 && y == 0)
return (cpackf(0, -y));
return (CMPLXF(0, -y));
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (cpackf(pio2_hi - (x - pio2_lo), -y));
return (CMPLXF(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
@ -253,7 +253,7 @@ cacosf(float complex z)
}
if (sy == 0)
ry = -ry;
return (cpackf(rx, ry));
return (CMPLXF(rx, ry));
}
float complex
@ -266,12 +266,12 @@ cacoshf(float complex z)
rx = crealf(w);
ry = cimagf(w);
if (isnan(rx) && isnan(ry))
return (cpackf(ry, rx));
return (CMPLXF(ry, rx));
if (isnan(rx))
return (cpackf(fabsf(ry), rx));
return (CMPLXF(fabsf(ry), rx));
if (isnan(ry))
return (cpackf(ry, ry));
return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
return (CMPLXF(ry, ry));
return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z))));
}
static float complex
@ -291,13 +291,13 @@ clog_for_large_values(float complex z)
}
if (ax > FLT_MAX / 2)
return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1,
atan2f(y, x)));
if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
}
static inline float
@ -346,22 +346,22 @@ catanhf(float complex z)
ay = fabsf(y);
if (y == 0 && ax <= 1)
return (cpackf(atanhf(x), y));
return (CMPLXF(atanhf(x), y));
if (x == 0)
return (cpackf(x, atanf(y)));
return (CMPLXF(x, atanf(y)));
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(copysignf(0, x), y + y));
return (CMPLXF(copysignf(0, x), y + y));
if (isinf(y))
return (cpackf(copysignf(0, x),
return (CMPLXF(copysignf(0, x),
copysignf(pio2_hi + pio2_lo, y)));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
return (cpackf(real_part_reciprocal(x, y),
return (CMPLXF(real_part_reciprocal(x, y),
copysignf(pio2_hi + pio2_lo, y)));
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@ -381,13 +381,13 @@ catanhf(float complex z)
else
ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
}
float complex
catanf(float complex z)
{
float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
return (CMPLXF(cimagf(w), crealf(w)));
}

30
lib/msun/src/e_j0.c
View File

@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static double pzero(double), qzero(double);
static __inline double pzero(double), qzero(double);
static const volatile double vone = 1, vzero = 0;
static const double
huge = 1e300,
@ -115,7 +117,7 @@ __ieee754_j0(double x)
if(ix<0x3f200000) { /* |x| < 2**-13 */
if(huge+x>one) { /* raise inexact if x != 0 */
if(ix<0x3e400000) return one; /* |x|<2**-27 */
else return one - 0.25*x*x;
else return one - x*x/4;
}
}
z = x*x;
@ -150,10 +152,16 @@ __ieee754_y0(double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/*
* y0(NaN) = NaN.
* y0(Inf) = 0.
* y0(-Inf) = NaN and raise invalid exception.
*/
if(ix>=0x7ff00000) return vone/(x+x*x);
/* y0(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y0(x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -268,7 +276,8 @@ static const double pS2[5] = {
1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
};
static double pzero(double x)
static __inline double
pzero(double x)
{
const double *p,*q;
double z,r,s;
@ -278,7 +287,7 @@ static const double pS2[5] = {
if(ix>=0x40200000) {p = pR8; q= pS8;}
else if(ix>=0x40122E8B){p = pR5; q= pS5;}
else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
else if(ix>=0x40000000){p = pR2; q= pS2;}
else {p = pR2; q= pS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -363,7 +372,8 @@ static const double qS2[6] = {
-5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
};
static double qzero(double x)
static __inline double
qzero(double x)
{
const double *p,*q;
double s,r,z;
@ -373,7 +383,7 @@ static const double qS2[6] = {
if(ix>=0x40200000) {p = qR8; q= qS8;}
else if(ix>=0x40122E8B){p = qR5; q= qS5;}
else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
else if(ix>=0x40000000){p = qR2; q= qS2;}
else {p = qR2; q= qS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

45
lib/msun/src/e_j0f.c
View File

@ -16,10 +16,16 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_j0.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static float pzerof(float), qzerof(float);
static __inline float pzerof(float), qzerof(float);
static const volatile float vone = 1, vzero = 0;
static const float
huge = 1e30,
@ -62,17 +68,17 @@ __ieee754_j0f(float x)
* 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>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
}
return z;
}
if(ix<0x39000000) { /* |x| < 2**-13 */
if(ix<0x3b000000) { /* |x| < 2**-9 */
if(huge+x>one) { /* raise inexact if x != 0 */
if(ix<0x32000000) return one; /* |x|<2**-27 */
else return one - (float)0.25*x*x;
if(ix<0x39800000) return one; /* |x|<2**-12 */
else return one - x*x/4;
}
}
z = x*x;
@ -107,10 +113,9 @@ __ieee754_y0f(float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x);
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix>=0x7f800000) return vone/(x+x*x);
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -136,14 +141,14 @@ __ieee754_y0f(float x)
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x32000000) { /* x < 2**-27 */
if(ix<=0x39000000) { /* x < 2**-13 */
return(u00 + tpi*__ieee754_logf(x));
}
z = x*x;
@ -224,7 +229,8 @@ static const float pS2[5] = {
1.4657617569e+01, /* 0x416a859a */
};
static float pzerof(float x)
static __inline float
pzerof(float x)
{
const float *p,*q;
float z,r,s;
@ -232,9 +238,9 @@ static const float pS2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pR8; q= pS8;}
else if(ix>=0x40f71c58){p = pR5; q= pS5;}
else if(ix>=0x4036db68){p = pR3; q= pS3;}
else if(ix>=0x40000000){p = pR2; q= pS2;}
else if(ix>=0x409173eb){p = pR5; q= pS5;}
else if(ix>=0x4036d917){p = pR3; q= pS3;}
else {p = pR2; q= pS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -319,7 +325,8 @@ static const float qS2[6] = {
-5.3109550476e+00, /* 0xc0a9f358 */
};
static float qzerof(float x)
static __inline float
qzerof(float x)
{
const float *p,*q;
float s,r,z;
@ -327,9 +334,9 @@ static const float qS2[6] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = qR8; q= qS8;}
else if(ix>=0x40f71c58){p = qR5; q= qS5;}
else if(ix>=0x4036db68){p = qR3; q= qS3;}
else if(ix>=0x40000000){p = qR2; q= qS2;}
else if(ix>=0x409173eb){p = qR5; q= qS5;}
else if(ix>=0x4036d917){p = qR3; q= qS3;}
else {p = qR2; q= qS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

28
lib/msun/src/e_j1.c
View File

@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static double pone(double), qone(double);
static __inline double pone(double), qone(double);
static const volatile double vone = 1, vzero = 0;
static const double
huge = 1e300,
@ -147,10 +149,16 @@ __ieee754_y1(double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/*
* y1(NaN) = NaN.
* y1(Inf) = 0.
* y1(-Inf) = NaN and raise invalid exception.
*/
if(ix>=0x7ff00000) return vone/(x+x*x);
/* y1(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y1(x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
@ -262,7 +270,8 @@ static const double ps2[5] = {
8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
};
static double pone(double x)
static __inline double
pone(double x)
{
const double *p,*q;
double z,r,s;
@ -272,7 +281,7 @@ static const double ps2[5] = {
if(ix>=0x40200000) {p = pr8; q= ps8;}
else if(ix>=0x40122E8B){p = pr5; q= ps5;}
else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
else if(ix>=0x40000000){p = pr2; q= ps2;}
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -358,7 +367,8 @@ static const double qs2[6] = {
-4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
};
static double qone(double x)
static __inline double
qone(double x)
{
const double *p,*q;
double s,r,z;
@ -368,7 +378,7 @@ static const double qs2[6] = {
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40122E8B){p = qr5; q= qs5;}
else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
else if(ix>=0x40000000){p = qr2; q= qs2;}
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

43
lib/msun/src/e_j1f.c
View File

@ -16,10 +16,16 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_j1.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static float ponef(float), qonef(float);
static __inline float ponef(float), qonef(float);
static const volatile float vone = 1, vzero = 0;
static const float
huge = 1e30,
@ -63,7 +69,7 @@ __ieee754_j1f(float x)
* 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>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
else {
u = ponef(y); v = qonef(y);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
@ -71,7 +77,7 @@ __ieee754_j1f(float x)
if(hx<0) return -z;
else return z;
}
if(ix<0x32000000) { /* |x|<2**-27 */
if(ix<0x39000000) { /* |x|<2**-13 */
if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
}
z = x*x;
@ -104,10 +110,9 @@ __ieee754_y1f(float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x);
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix>=0x7f800000) return vone/(x+x*x);
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
@ -129,14 +134,14 @@ __ieee754_y1f(float x)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = ponef(x); v = qonef(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x24800000) { /* x < 2**-54 */
if(ix<=0x33000000) { /* x < 2**-25 */
return(-tpi/x);
}
z = x*x;
@ -219,7 +224,8 @@ static const float ps2[5] = {
8.3646392822e+00, /* 0x4105d590 */
};
static float ponef(float x)
static __inline float
ponef(float x)
{
const float *p,*q;
float z,r,s;
@ -227,9 +233,9 @@ static const float ps2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pr8; q= ps8;}
else if(ix>=0x40f71c58){p = pr5; q= ps5;}
else if(ix>=0x4036db68){p = pr3; q= ps3;}
else if(ix>=0x40000000){p = pr2; q= ps2;}
else if(ix>=0x409173eb){p = pr5; q= ps5;}
else if(ix>=0x4036d917){p = pr3; q= ps3;}
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -315,17 +321,18 @@ static const float qs2[6] = {
-4.9594988823e+00, /* 0xc09eb437 */
};
static float qonef(float x)
static __inline float
qonef(float x)
{
const float *p,*q;
float s,r,z;
int32_t ix;
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40f71c58){p = qr5; q= qs5;}
else if(ix>=0x4036db68){p = qr3; q= qs3;}
else if(ix>=0x40000000){p = qr2; q= qs2;}
if(ix>=0x41000000) {p = qr8; q= qs8;}
else if(ix>=0x409173eb){p = qr5; q= qs5;}
else if(ix>=0x4036d917){p = qr3; q= qs3;}
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

10
lib/msun/src/e_jn.c
View File

@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static const volatile double vone = 1, vzero = 0;
static const double
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
@ -220,10 +222,12 @@ __ieee754_yn(int n, double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
/* yn(n,NaN) = NaN */
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/* yn(n,+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* yn(n,x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
sign = 1;
if(n<0){
n = -n;

11
lib/msun/src/e_jnf.c
View File

@ -16,9 +16,15 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_jn.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static const volatile float vone = 1, vzero = 0;
static const float
two = 2.0000000000e+00, /* 0x40000000 */
one = 1.0000000000e+00; /* 0x3F800000 */
@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
if(ix>0x7f800000) return x+x;
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
sign = 1;
if(n<0){
n = -n;

6
lib/msun/src/e_lgamma.c
View File

@ -21,6 +21,8 @@ __FBSDID("$FreeBSD$");
* Method: call __ieee754_lgamma_r
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -31,3 +33,7 @@ __ieee754_lgamma(double x)
{
return __ieee754_lgamma_r(x,&signgam);
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(lgamma, lgammal);
#endif

88
lib/msun/src/e_lgamma_r.c
View File

@ -1,4 +1,3 @@
/* @(#)e_lgamma_r.c 1.3 95/01/18 */
/*
* ====================================================
@ -6,22 +5,21 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_lgamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
@ -59,20 +57,20 @@ __FBSDID("$FreeBSD$");
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
* where
* where
* |w - f(z)| < 2**-58.74
*
*
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
* Note: one should avoid compute pi*(-x) directly in the
* Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
*
*
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1) = lgamma(2) = 0
@ -80,9 +78,10 @@ __FBSDID("$FreeBSD$");
* lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
* lgamma(inf) = inf
* lgamma(-inf) = inf (bug for bug compatible with C99!?)
*
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -187,9 +186,9 @@ sin_pi(double x)
switch (n) {
case 0: y = __kernel_sin(pi*y,zero,0); break;
case 1:
case 1:
case 2: y = __kernel_cos(pi*(0.5-y),zero); break;
case 3:
case 3:
case 4: y = __kernel_sin(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cos(pi*(y-1.5),zero); break;
@ -202,24 +201,28 @@ sin_pi(double x)
double
__ieee754_lgamma_r(double x, int *signgamp)
{
double t,y,z,nadj,p,p1,p2,p3,q,r,w;
double nadj,p,p1,p2,p3,q,r,t,w,y,z;
int32_t hx;
int i,ix,lx;
EXTRACT_WORDS(hx,lx,x);
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
/* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x*x;
if((ix|lx)==0) return one/vzero;
if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx<0) {
*signgamp = -1;
return -__ieee754_log(-x);
} else return -__ieee754_log(x);
/* purge +-0 and tiny arguments */
*signgamp = 1-2*((uint32_t)hx>>31);
if(ix<0x3c700000) { /* |x|<2**-56, return -log(|x|) */
if((ix|lx)==0)
return one/vzero;
return -__ieee754_log(fabs(x));
}
/* purge negative integers and start evaluation for other x < 0 */
if(hx<0) {
*signgamp = 1;
if(ix>=0x43300000) /* |x|>=2**52, must be -integer */
return one/vzero;
t = sin_pi(x);
@ -229,7 +232,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
x = -x;
}
/* purge off 1 and 2 */
/* purge 1 and 2 */
if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
/* for x < 2.0 */
else if(ix<0x40000000) {
@ -250,7 +253,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
p = y*p1+p2;
r += (p-0.5*y); break;
r += p-y/2; break;
case 1:
z = y*y;
w = z*y;
@ -258,38 +261,43 @@ __ieee754_lgamma_r(double x, int *signgamp)
p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
p = z*p1-(tt-w*(p2+y*p3));
r += (tf + p); break;
case 2:
r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += (-0.5*y + p1/p2);
r += p1/p2-y/2;
}
}
else if(ix<0x40200000) { /* x < 8.0 */
i = (int)x;
y = x-(double)i;
/* x < 8.0 */
else if(ix<0x40200000) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = half*y+p/q;
r = y/2+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+6.0); /* FALLTHRU */
case 6: z *= (y+5.0); /* FALLTHRU */
case 5: z *= (y+4.0); /* FALLTHRU */
case 4: z *= (y+3.0); /* FALLTHRU */
case 3: z *= (y+2.0); /* FALLTHRU */
case 7: z *= (y+6); /* FALLTHRU */
case 6: z *= (y+5); /* FALLTHRU */
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_log(z); break;
}
/* 8.0 <= x < 2**58 */
} else if (ix < 0x43900000) {
/* 8.0 <= x < 2**56 */
} else if (ix < 0x43700000) {
t = __ieee754_log(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
} else
/* 2**58 <= x <= inf */
} else
/* 2**56 <= x <= inf */
r = x*(__ieee754_log(x)-one);
if(hx<0) r = nadj - r;
return r;
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(lgamma_r, lgammal_r);
#endif

200
lib/msun/src/e_lgammaf_r.c
View File

@ -1,5 +1,6 @@
/* e_lgammaf_r.c -- float version of e_lgamma_r.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Conversion to float fixed By Steven G. Kargl.
*/
/*
@ -22,72 +23,63 @@ __FBSDID("$FreeBSD$");
static const volatile float vzero = 0;
static const float
zero= 0.0000000000e+00,
half= 5.0000000000e-01, /* 0x3f000000 */
one = 1.0000000000e+00, /* 0x3f800000 */
zero= 0,
half= 0.5,
one = 1,
pi = 3.1415927410e+00, /* 0x40490fdb */
a0 = 7.7215664089e-02, /* 0x3d9e233f */
a1 = 3.2246702909e-01, /* 0x3ea51a66 */
a2 = 6.7352302372e-02, /* 0x3d89f001 */
a3 = 2.0580807701e-02, /* 0x3ca89915 */
a4 = 7.3855509982e-03, /* 0x3bf2027e */
a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */
a6 = 1.1927076848e-03, /* 0x3a9c54a1 */
a7 = 5.1006977446e-04, /* 0x3a05b634 */
a8 = 2.2086278477e-04, /* 0x39679767 */
a9 = 1.0801156895e-04, /* 0x38e28445 */
a10 = 2.5214456400e-05, /* 0x37d383a2 */
a11 = 4.4864096708e-05, /* 0x383c2c75 */
tc = 1.4616321325e+00, /* 0x3fbb16c3 */
tf = -1.2148628384e-01, /* 0xbdf8cdcd */
/* tt = -(tail of tf) */
tt = 6.6971006518e-09, /* 0x31e61c52 */
t0 = 4.8383611441e-01, /* 0x3ef7b95e */
t1 = -1.4758771658e-01, /* 0xbe17213c */
t2 = 6.4624942839e-02, /* 0x3d845a15 */
t3 = -3.2788541168e-02, /* 0xbd064d47 */
t4 = 1.7970675603e-02, /* 0x3c93373d */
t5 = -1.0314224288e-02, /* 0xbc28fcfe */
t6 = 6.1005386524e-03, /* 0x3bc7e707 */
t7 = -3.6845202558e-03, /* 0xbb7177fe */
t8 = 2.2596477065e-03, /* 0x3b141699 */
t9 = -1.4034647029e-03, /* 0xbab7f476 */
t10 = 8.8108185446e-04, /* 0x3a66f867 */
t11 = -5.3859531181e-04, /* 0xba0d3085 */
t12 = 3.1563205994e-04, /* 0x39a57b6b */
t13 = -3.1275415677e-04, /* 0xb9a3f927 */
t14 = 3.3552918467e-04, /* 0x39afe9f7 */
u0 = -7.7215664089e-02, /* 0xbd9e233f */
u1 = 6.3282704353e-01, /* 0x3f2200f4 */
u2 = 1.4549225569e+00, /* 0x3fba3ae7 */
u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */
u4 = 2.2896373272e-01, /* 0x3e6a7578 */
u5 = 1.3381091878e-02, /* 0x3c5b3c5e */
v1 = 2.4559779167e+00, /* 0x401d2ebe */
v2 = 2.1284897327e+00, /* 0x4008392d */
v3 = 7.6928514242e-01, /* 0x3f44efdf */
v4 = 1.0422264785e-01, /* 0x3dd572af */
v5 = 3.2170924824e-03, /* 0x3b52d5db */
s0 = -7.7215664089e-02, /* 0xbd9e233f */
s1 = 2.1498242021e-01, /* 0x3e5c245a */
s2 = 3.2577878237e-01, /* 0x3ea6cc7a */
s3 = 1.4635047317e-01, /* 0x3e15dce6 */
s4 = 2.6642270386e-02, /* 0x3cda40e4 */
s5 = 1.8402845599e-03, /* 0x3af135b4 */
s6 = 3.1947532989e-05, /* 0x3805ff67 */
r1 = 1.3920053244e+00, /* 0x3fb22d3b */
r2 = 7.2193557024e-01, /* 0x3f38d0c5 */
r3 = 1.7193385959e-01, /* 0x3e300f6e */
r4 = 1.8645919859e-02, /* 0x3c98bf54 */
r5 = 7.7794247773e-04, /* 0x3a4beed6 */
r6 = 7.3266842264e-06, /* 0x36f5d7bd */
w0 = 4.1893854737e-01, /* 0x3ed67f1d */
w1 = 8.3333335817e-02, /* 0x3daaaaab */
w2 = -2.7777778450e-03, /* 0xbb360b61 */
w3 = 7.9365057172e-04, /* 0x3a500cfd */
w4 = -5.9518753551e-04, /* 0xba1c065c */
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
/*
* Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]:
* |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4
*/
a0 = 7.72156641e-02, /* 0x3d9e233f */
a1 = 3.22467119e-01, /* 0x3ea51a69 */
a2 = 6.73484802e-02, /* 0x3d89ee00 */
a3 = 2.06395667e-02, /* 0x3ca9144f */
a4 = 6.98275631e-03, /* 0x3be4cf9b */
a5 = 4.11768444e-03, /* 0x3b86eda4 */
/*
* Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]:
* |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8.
*/
tc = 1.46163213e+00, /* 0x3fbb16c3 */
tf = -1.21486291e-01, /* 0xbdf8cdce */
t0 = -2.94064460e-11, /* 0xae0154b7 */
t1 = -2.35939837e-08, /* 0xb2caabb8 */
t2 = 4.83836412e-01, /* 0x3ef7b968 */
t3 = -1.47586212e-01, /* 0xbe1720d7 */
t4 = 6.46013096e-02, /* 0x3d844db1 */
t5 = -3.28450352e-02, /* 0xbd068884 */
t6 = 1.86483748e-02, /* 0x3c98c47a */
t7 = -9.89206228e-03, /* 0xbc221251 */
/*
* Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]:
* |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2
*/
u0 = -7.72156641e-02, /* 0xbd9e233f */
u1 = 7.36789703e-01, /* 0x3f3c9e40 */
u2 = 4.95649040e-01, /* 0x3efdc5b6 */
v1 = 1.10958421e+00, /* 0x3f8e06db */
v2 = 2.10598111e-01, /* 0x3e57a708 */
v3 = -1.02995494e-02, /* 0xbc28bf71 */
/*
* Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]:
* |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0
* with y = x - 2.
*/
s0 = -7.72156641e-02, /* 0xbd9e233f */
s1 = 2.69987404e-01, /* 0x3e8a3bca */
s2 = 1.42851010e-01, /* 0x3e124789 */
s3 = 1.19389519e-02, /* 0x3c439b98 */
r1 = 6.79650068e-01, /* 0x3f2dfd8c */
r2 = 1.16058730e-01, /* 0x3dedb033 */
r3 = 3.75673687e-03, /* 0x3b763396 */
/*
* Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]:
* |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6.
*/
w0 = 4.18938547e-01, /* 0x3ed67f1d */
w1 = 8.33332464e-02, /* 0x3daaaa9f */
w2 = -2.76129087e-03; /* 0xbb34f6c6 */
static float
sin_pif(float x)
@ -130,25 +122,29 @@ sin_pif(float x)
float
__ieee754_lgammaf_r(float x, int *signgamp)
{
float t,y,z,nadj,p,p1,p2,p3,q,r,w;
float nadj,p,p1,p2,p3,q,r,t,w,y,z;
int32_t hx;
int i,ix;
GET_FLOAT_WORD(hx,x);
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
/* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7f800000) return x*x;
if(ix==0) return one/vzero;
if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */
if(hx<0) {
*signgamp = -1;
return -__ieee754_logf(-x);
} else return -__ieee754_logf(x);
/* purge +-0 and tiny arguments */
*signgamp = 1-2*((uint32_t)hx>>31);
if(ix<0x32000000) { /* |x|<2**-27, return -log(|x|) */
if(ix==0)
return one/vzero;
return -__ieee754_logf(fabsf(x));
}
/* purge negative integers and start evaluation for other x < 0 */
if(hx<0) {
if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
*signgamp = 1;
if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
return one/vzero;
t = sin_pif(x);
if(t==zero) return one/vzero; /* -integer */
@ -157,7 +153,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
x = -x;
}
/* purge off 1 and 2 */
/* purge 1 and 2 */
if (ix==0x3f800000||ix==0x40000000) r = 0;
/* for x < 2.0 */
else if(ix<0x40000000) {
@ -168,55 +164,51 @@ __ieee754_lgammaf_r(float x, int *signgamp)
else {y = x; i=2;}
} else {
r = zero;
if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
if(ix>=0x3fdda618) {y=2-x;i=0;} /* [1.7316,2] */
else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
else {y=x-one;i=2;}
}
switch(i) {
case 0:
z = y*y;
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
p1 = a0+z*(a2+z*a4);
p2 = z*(a1+z*(a3+z*a5));
p = y*p1+p2;
r += (p-(float)0.5*y); break;
r += p-y/2; break;
case 1:
z = y*y;
w = z*y;
p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
p = z*p1-(tt-w*(p2+y*p3));
r += (tf + p); break;
p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7)))));
r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += (-(float)0.5*y + p1/p2);
p1 = y*(u0+y*(u1+y*u2));
p2 = one+y*(v1+y*(v2+y*v3));
r += p1/p2-y/2;
}
}
else if(ix<0x41000000) { /* x < 8.0 */
i = (int)x;
y = x-(float)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = half*y+p/q;
/* x < 8.0 */
else if(ix<0x41000000) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*s3)));
q = one+y*(r1+y*(r2+y*r3));
r = y/2+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+(float)6.0); /* FALLTHRU */
case 6: z *= (y+(float)5.0); /* FALLTHRU */
case 5: z *= (y+(float)4.0); /* FALLTHRU */
case 4: z *= (y+(float)3.0); /* FALLTHRU */
case 3: z *= (y+(float)2.0); /* FALLTHRU */
case 7: z *= (y+6); /* FALLTHRU */
case 6: z *= (y+5); /* FALLTHRU */
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_logf(z); break;
}
/* 8.0 <= x < 2**58 */
} else if (ix < 0x5c800000) {
/* 8.0 <= x < 2**27 */
} else if (ix < 0x4d000000) {
t = __ieee754_logf(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
w = w0+z*(w1+y*w2);
r = (x-half)*(t-one)+w;
} else
/* 2**58 <= x <= inf */
/* 2**27 <= x <= inf */
r = x*(__ieee754_logf(x)-one);
if(hx<0) r = nadj - r;
return r;

25
lib/msun/src/e_lgammal.c Normal file
View File

@ -0,0 +1,25 @@
/* @(#)e_lgamma.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
extern int signgam;
long double
lgammal(long double x)
{
return lgammal_r(x,&signgam);
}

1
lib/msun/src/imprecise.c
View File

@ -60,5 +60,4 @@ DECLARE_WEAK(powl);
long double imprecise_ ## f ## l(long double v) { return f(v); }\
DECLARE_WEAK(f ## l)
DECLARE_IMPRECISE(lgamma);
DECLARE_IMPRECISE(tgamma);

2
lib/msun/src/k_exp.c
View File

@ -103,6 +103,6 @@ __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
return (cpack(cos(y) * exp_x * scale1 * scale2,
return (CMPLX(cos(y) * exp_x * scale1 * scale2,
sin(y) * exp_x * scale1 * scale2));
}

2
lib/msun/src/k_expf.c
View File

@ -82,6 +82,6 @@ __ldexp_cexpf(float complex z, int expt)
half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
return (cpackf(cosf(y) * exp_x * scale1 * scale2,
return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,
sinf(y) * exp_x * scale1 * scale2));
}

6
lib/msun/src/math.h
View File

@ -500,8 +500,12 @@ long double tanhl(long double);
long double tanl(long double);
long double tgammal(long double);
long double truncl(long double);
#endif /* __ISO_C_VISIBLE >= 1999 */
#if __BSD_VISIBLE
long double lgammal_r(long double, int *);
#endif
__END_DECLS
#endif /* !_MATH_H_ */

18
lib/msun/src/math_private.h
View File

@ -454,9 +454,15 @@ typedef union {
* (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
* In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
* to -0.0+I*0.0.
*
* The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
* to construct complex values. Compilers that conform to the C99
* standard require the following functions to avoid the above issues.
*/
#ifndef CMPLXF
static __inline float complex
cpackf(float x, float y)
CMPLXF(float x, float y)
{
float_complex z;
@ -464,9 +470,11 @@ cpackf(float x, float y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#ifndef CMPLX
static __inline double complex
cpack(double x, double y)
CMPLX(double x, double y)
{
double_complex z;
@ -474,9 +482,11 @@ cpack(double x, double y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#ifndef CMPLXL
static __inline long double complex
cpackl(long double x, long double y)
CMPLXL(long double x, long double y)
{
long_double_complex z;
@ -484,6 +494,8 @@ cpackl(long double x, long double y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#endif /* _COMPLEX_H */
#ifdef __GNUCLIKE_ASM

59
lib/msun/src/s_ccosh.c
View File

@ -32,6 +32,8 @@
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
* The sign of the result for some exceptional values is unspecified but
* must satisfy both cosh(conj(z)) == conj(cosh(z)) and cosh(-z) == cosh(z).
*/
#include <sys/cdefs.h>
@ -62,49 +64,48 @@ ccosh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return (cpack(cosh(x), x * y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
return (CMPLX(cosh(x), x * y));
if (ix < 0x40360000) /* |x| < 22: normal case */
return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return (cpack(h * cos(y), copysign(h, x) * sin(y)));
return (CMPLX(h * cos(y), copysign(h, x) * sin(y)));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return (cpack(creal(z), cimag(z) * copysign(1, x)));
z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return (CMPLX(creal(z), cimag(z) * copysign(1, x)));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return (cpack(h * h * cos(y), h * sin(y)));
return (CMPLX(h * h * cos(y), h * sin(y)));
}
}
/*
* cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
* cosh(+-0 +- I Inf) = dNaN + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument. Raise the invalid floating-point
* exception.
*
* cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
* cosh(+-0 +- I NaN) = d(NaN) + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument.
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(y - y, copysign(0, x * (y - y))));
if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, x * copysign(0, y)));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
*
* cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
* The sign of 0 in the result is unspecified.
* cosh(NaN +- I 0) = d(NaN) + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x * x, copysign(0, x) * y));
return (cpack(x * x, copysign(0, (x + x) * y)));
}
if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */
return (CMPLX(x * x, copysign(0, x) * y));
/*
* cosh(x +- I Inf) = dNaN + I dNaN.
@ -114,8 +115,8 @@ ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, x * (y - y)));
/*
* cosh(+-Inf + I NaN) = +Inf + I d(NaN).
@ -126,10 +127,10 @@ ccosh(double complex z)
*
* cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (ix == 0x7ff00000 && lx == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack((x * x) * cos(y), x * sin(y)));
return (CMPLX(INFINITY, x * (y - y)));
return (CMPLX(INFINITY * cos(y), x * sin(y)));
}
/*
@ -143,7 +144,7 @@ ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (cpack((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
}
double complex
@ -151,5 +152,5 @@ ccos(double complex z)
{
/* ccos(z) = ccosh(I * z) */
return (ccosh(cpack(-cimag(z), creal(z))));
return (ccosh(CMPLX(-cimag(z), creal(z))));
}

43
lib/msun/src/s_ccoshf.c
View File

@ -25,7 +25,7 @@
*/
/*
* Hyperbolic cosine of a complex argument. See s_ccosh.c for details.
* Float version of ccosh(). See s_ccosh.c for details.
*/
#include <sys/cdefs.h>
@ -55,50 +55,47 @@ ccoshf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(coshf(x), x * y));
if (ix < 0x41100000) /* small x: normal case */
return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
return (CMPLXF(coshf(x), x * y));
if (ix < 0x41100000) /* |x| < 9: normal case */
return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return (cpackf(h * cosf(y), copysignf(h, x) * sinf(y)));
h = expf(fabsf(x)) * 0.5F;
return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y)));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return (cpackf(crealf(z), cimagf(z) * copysignf(1, x)));
z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x)));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return (cpackf(h * h * cosf(y), h * sinf(y)));
return (CMPLXF(h * h * cosf(y), h * sinf(y)));
}
}
if (ix == 0 && iy >= 0x7f800000)
return (cpackf(y - y, copysignf(0, x * (y - y))));
if (ix == 0) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, x * copysignf(0, y)));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return (cpackf(x * x, copysignf(0, x) * y));
return (cpackf(x * x, copysignf(0, (x + x) * y)));
}
if (iy == 0) /* && ix >= 0x7f800000 */
return (CMPLXF(x * x, copysignf(0, x) * y));
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
if (ix < 0x7f800000) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, x * (y - y)));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (ix == 0x7f800000) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf((x * x) * cosf(y), x * sinf(y)));
return (CMPLXF(INFINITY, x * (y - y)));
return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
}
return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
}
float complex
ccosf(float complex z)
{
return (ccoshf(cpackf(-cimagf(z), crealf(z))));
return (ccoshf(CMPLXF(-cimagf(z), crealf(z))));
}

12
lib/msun/src/s_cexp.c
View File

@ -50,22 +50,22 @@ cexp(double complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if ((hy | ly) == 0)
return (cpack(exp(x), y));
return (CMPLX(exp(x), y));
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0)
return (cpack(cos(y), sin(y)));
return (CMPLX(cos(y), sin(y)));
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return (cpack(y - y, y - y));
return (CMPLX(y - y, y - y));
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return (cpack(0.0, 0.0));
return (CMPLX(0.0, 0.0));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return (cpack(x, y - y));
return (CMPLX(x, y - y));
}
}
@ -84,6 +84,6 @@ cexp(double complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = exp(x);
return (cpack(exp_x * cos(y), exp_x * sin(y)));
return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
}
}

12
lib/msun/src/s_cexpf.c
View File

@ -50,22 +50,22 @@ cexpf(float complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if (hy == 0)
return (cpackf(expf(x), y));
return (CMPLXF(expf(x), y));
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0)
return (cpackf(cosf(y), sinf(y)));
return (CMPLXF(cosf(y), sinf(y)));
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return (cpackf(y - y, y - y));
return (CMPLXF(y - y, y - y));
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return (cpackf(0.0, 0.0));
return (CMPLXF(0.0, 0.0));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return (cpackf(x, y - y));
return (CMPLXF(x, y - y));
}
}
@ -84,6 +84,6 @@ cexpf(float complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = expf(x);
return (cpackf(exp_x * cosf(y), exp_x * sinf(y)));
return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
}
}

2
lib/msun/src/s_conj.c
View File

@ -34,5 +34,5 @@ double complex
conj(double complex z)
{
return (cpack(creal(z), -cimag(z)));
return (CMPLX(creal(z), -cimag(z)));
}

2
lib/msun/src/s_conjf.c
View File

@ -34,5 +34,5 @@ float complex
conjf(float complex z)
{
return (cpackf(crealf(z), -cimagf(z)));
return (CMPLXF(crealf(z), -cimagf(z)));
}

2
lib/msun/src/s_conjl.c
View File

@ -34,5 +34,5 @@ long double complex
conjl(long double complex z)
{
return (cpackl(creall(z), -cimagl(z)));
return (CMPLXL(creall(z), -cimagl(z)));
}

2
lib/msun/src/s_cproj.c
View File

@ -39,7 +39,7 @@ cproj(double complex z)
if (!isinf(creal(z)) && !isinf(cimag(z)))
return (z);
else
return (cpack(INFINITY, copysign(0.0, cimag(z))));
return (CMPLX(INFINITY, copysign(0.0, cimag(z))));
}
#if LDBL_MANT_DIG == 53

2
lib/msun/src/s_cprojf.c
View File

@ -39,5 +39,5 @@ cprojf(float complex z)
if (!isinf(crealf(z)) && !isinf(cimagf(z)))
return (z);
else
return (cpackf(INFINITY, copysignf(0.0, cimagf(z))));
return (CMPLXF(INFINITY, copysignf(0.0, cimagf(z))));
}

2
lib/msun/src/s_cprojl.c
View File

@ -39,5 +39,5 @@ cprojl(long double complex z)
if (!isinf(creall(z)) && !isinf(cimagl(z)))
return (z);
else
return (cpackl(INFINITY, copysignl(0.0, cimagl(z))));
return (CMPLXL(INFINITY, copysignl(0.0, cimagl(z))));
}

73
lib/msun/src/s_csinh.c
View File

@ -32,6 +32,8 @@
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
* The sign of the result for some exceptional values is unspecified but
* must satisfy both sinh(conj(z)) == conj(sinh(z)) and sinh(-z) == -sinh(z).
*/
#include <sys/cdefs.h>
@ -62,48 +64,45 @@ csinh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return (cpack(sinh(x), y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(sinh(x) * cos(y), cosh(x) * sin(y)));
return (CMPLX(sinh(x), y));
if (ix < 0x40360000) /* |x| < 22: normal case */
return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y)));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return (cpack(copysign(h, x) * cos(y), h * sin(y)));
return (CMPLX(copysign(h, x) * cos(y), h * sin(y)));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return (cpack(creal(z) * copysign(1, x), cimag(z)));
z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return (CMPLX(creal(z) * copysign(1, x), cimag(z)));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return (cpack(h * cos(y), h * h * sin(y)));
return (CMPLX(h * cos(y), h * h * sin(y)));
}
}
/*
* sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
* sinh(+-0 +- I Inf) = +-0 + I dNaN.
* The sign of 0 in the result is unspecified. Choice = same sign
* as the argument. Raise the invalid floating-point exception.
*
* sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
* sinh(+-0 +- I NaN) = +-0 + I d(NaN).
* The sign of 0 in the result is unspecified. Choice = same sign
* as the argument.
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(copysign(0, x * (y - y)), y - y));
if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */
return (CMPLX(x, y - y));
/*
* sinh(+-Inf +- I 0) = +-Inf + I +-0.
*
* sinh(NaN +- I 0) = d(NaN) + I +-0.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x, y));
return (cpack(x, copysign(0, y)));
}
if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */
return (CMPLX(x + x, y));
/*
* sinh(x +- I Inf) = dNaN + I dNaN.
@ -113,45 +112,45 @@ csinh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, y - y));
/*
* sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
* The sign of Inf in the result is unspecified. Choice = normally
* the same as d(NaN).
* The sign of Inf in the result is unspecified. Choice = same sign
* as the argument.
*
* sinh(+-Inf +- I Inf) = +Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = always +.
* Raise the invalid floating-point exception.
* sinh(+-Inf +- I Inf) = +-Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = same sign
* as the argument. Raise the invalid floating-point exception.
*
* sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (ix == 0x7ff00000 && lx == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack(x * cos(y), INFINITY * sin(y)));
return (CMPLX(x, y - y));
return (CMPLX(x * cos(y), INFINITY * sin(y)));
}
/*
* sinh(NaN + I NaN) = d(NaN) + I d(NaN).
* sinh(NaN1 + I NaN2) = d(NaN1, NaN2) + I d(NaN1, NaN2).
*
* sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
* sinh(NaN +- I Inf) = d(NaN, dNaN) + I d(NaN, dNaN).
* Optionally raises the invalid floating-point exception.
* Choice = raise.
*
* sinh(NaN + I y) = d(NaN) + I d(NaN).
* sinh(NaN + I y) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (cpack((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLX((x + x) * (y - y), (x * x) * (y - y)));
}
double complex
csin(double complex z)
{
/* csin(z) = -I * csinh(I * z) */
z = csinh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
/* csin(z) = -I * csinh(I * z) = I * conj(csinh(I * conj(z))). */
z = csinh(CMPLX(cimag(z), creal(z)));
return (CMPLX(cimag(z), creal(z)));
}

45
lib/msun/src/s_csinhf.c
View File

@ -25,7 +25,7 @@
*/
/*
* Hyperbolic sine of a complex argument z. See s_csinh.c for details.
* Float version of csinh(). See s_csinh.c for details.
*/
#include <sys/cdefs.h>
@ -55,51 +55,48 @@ csinhf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(sinhf(x), y));
if (ix < 0x41100000) /* small x: normal case */
return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
return (CMPLXF(sinhf(x), y));
if (ix < 0x41100000) /* |x| < 9: normal case */
return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y)));
h = expf(fabsf(x)) * 0.5F;
return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y)));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return (cpackf(crealf(z) * copysignf(1, x), cimagf(z)));
z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return (CMPLXF(crealf(z) * copysignf(1, x), cimagf(z)));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return (cpackf(h * cosf(y), h * h * sinf(y)));
return (CMPLXF(h * cosf(y), h * h * sinf(y)));
}
}
if (ix == 0 && iy >= 0x7f800000)
return (cpackf(copysignf(0, x * (y - y)), y - y));
if (ix == 0) /* && iy >= 0x7f800000 */
return (CMPLXF(x, y - y));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return (cpackf(x, y));
return (cpackf(x, copysignf(0, y)));
}
if (iy == 0) /* && ix >= 0x7f800000 */
return (CMPLXF(x + x, y));
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
if (ix < 0x7f800000) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (ix == 0x7f800000) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf(x * cosf(y), INFINITY * sinf(y)));
return (CMPLXF(x, y - y));
return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
}
return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLXF((x + x) * (y - y), (x * x) * (y - y)));
}
float complex
csinf(float complex z)
{
z = csinhf(cpackf(-cimagf(z), crealf(z)));
return (cpackf(cimagf(z), -crealf(z)));
z = csinhf(CMPLXF(cimagf(z), crealf(z)));
return (CMPLXF(cimagf(z), crealf(z)));
}

14
lib/msun/src/s_csqrt.c
View File

@ -58,12 +58,12 @@ csqrt(double complex z)
/* Handle special cases. */
if (z == 0)
return (cpack(0, b));
return (CMPLX(0, b));
if (isinf(b))
return (cpack(INFINITY, b));
return (CMPLX(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpack(a, t)); /* return NaN + NaN i */
return (CMPLX(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -73,9 +73,9 @@ csqrt(double complex z)
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpack(fabs(b - b), copysign(a, b)));
return (CMPLX(fabs(b - b), copysign(a, b)));
else
return (cpack(a, copysign(b - b, b)));
return (CMPLX(a, copysign(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -94,10 +94,10 @@ csqrt(double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
result = cpack(t, b / (2 * t));
result = CMPLX(t, b / (2 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
result = cpack(fabs(b) / (2 * t), copysign(t, b));
result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
}
/* Rescale. */

14
lib/msun/src/s_csqrtf.c
View File

@ -49,12 +49,12 @@ csqrtf(float complex z)
/* Handle special cases. */
if (z == 0)
return (cpackf(0, b));
return (CMPLXF(0, b));
if (isinf(b))
return (cpackf(INFINITY, b));
return (CMPLXF(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpackf(a, t)); /* return NaN + NaN i */
return (CMPLXF(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -64,9 +64,9 @@ csqrtf(float complex z)
* csqrtf(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpackf(fabsf(b - b), copysignf(a, b)));
return (CMPLXF(fabsf(b - b), copysignf(a, b)));
else
return (cpackf(a, copysignf(b - b, b)));
return (CMPLXF(a, copysignf(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -80,9 +80,9 @@ csqrtf(float complex z)
*/
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
return (cpackf(t, b / (2.0 * t)));
return (CMPLXF(t, b / (2.0 * t)));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b)));
return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)));
}
}

14
lib/msun/src/s_csqrtl.c
View File

@ -58,12 +58,12 @@ csqrtl(long double complex z)
/* Handle special cases. */
if (z == 0)
return (cpackl(0, b));
return (CMPLXL(0, b));
if (isinf(b))
return (cpackl(INFINITY, b));
return (CMPLXL(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpackl(a, t)); /* return NaN + NaN i */
return (CMPLXL(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -73,9 +73,9 @@ csqrtl(long double complex z)
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpackl(fabsl(b - b), copysignl(a, b)));
return (CMPLXL(fabsl(b - b), copysignl(a, b)));
else
return (cpackl(a, copysignl(b - b, b)));
return (CMPLXL(a, copysignl(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -94,10 +94,10 @@ csqrtl(long double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrtl((a + hypotl(a, b)) * 0.5);
result = cpackl(t, b / (2 * t));
result = CMPLXL(t, b / (2 * t));
} else {
t = sqrtl((-a + hypotl(a, b)) * 0.5);
result = cpackl(fabsl(b) / (2 * t), copysignl(t, b));
result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b));
}
/* Rescale. */

43
lib/msun/src/s_ctanh.c
View File

@ -25,7 +25,7 @@
*/
/*
* Hyperbolic tangent of a complex argument z = x + i y.
* Hyperbolic tangent of a complex argument z = x + I y.
*
* The algorithm is from:
*
@ -44,15 +44,15 @@
*
* tanh(z) = sinh(z) / cosh(z)
*
* sinh(x) cos(y) + i cosh(x) sin(y)
* sinh(x) cos(y) + I cosh(x) sin(y)
* = ---------------------------------
* cosh(x) cos(y) + i sinh(x) sin(y)
* cosh(x) cos(y) + I sinh(x) sin(y)
*
* cosh(x) sinh(x) / cos^2(y) + i tan(y)
* cosh(x) sinh(x) / cos^2(y) + I tan(y)
* = -------------------------------------
* 1 + sinh^2(x) / cos^2(y)
*
* beta rho s + i t
* beta rho s + I t
* = ----------------
* 1 + beta s^2
*
@ -85,16 +85,16 @@ ctanh(double complex z)
ix = hx & 0x7fffffff;
/*
* ctanh(NaN + i 0) = NaN + i 0
* ctanh(NaN +- I 0) = d(NaN) +- I 0
*
* ctanh(NaN + i y) = NaN + i NaN for y != 0
* ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0
*
* The imaginary part has the sign of x*sin(2*y), but there's no
* special effort to get this right.
*
* ctanh(+-Inf +- i Inf) = +-1 +- 0
* ctanh(+-Inf +- I Inf) = +-1 +- I 0
*
* ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
* ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite
*
* The imaginary part of the sign is unspecified. This special
* case is only needed to avoid a spurious invalid exception when
@ -102,26 +102,27 @@ ctanh(double complex z)
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
return (cpack(x, (y == 0 ? y : x * y)));
return (CMPLX((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
}
/*
* ctanh(x + i NAN) = NaN + i NaN
* ctanh(x +- i Inf) = NaN + i NaN
* ctanh(x + I NaN) = d(NaN) + I d(NaN)
* ctanh(x +- I Inf) = dNaN + I dNaN
*/
if (!isfinite(y))
return (cpack(y - y, y - y));
return (CMPLX(y - y, y - y));
/*
* ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
* ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the
* approximation sinh^2(huge) ~= exp(2*huge) / 4.
* We use a modified formula to avoid spurious overflow.
*/
if (ix >= 0x40360000) { /* x >= 22 */
if (ix >= 0x40360000) { /* |x| >= 22 */
double exp_mx = exp(-fabs(x));
return (cpack(copysign(1, x),
return (CMPLX(copysign(1, x),
4 * sin(y) * cos(y) * exp_mx * exp_mx));
}
@ -131,14 +132,14 @@ ctanh(double complex z)
s = sinh(x);
rho = sqrt(1 + s * s); /* = cosh(x) */
denom = 1 + beta * s * s;
return (cpack((beta * rho * s) / denom, t / denom));
return (CMPLX((beta * rho * s) / denom, t / denom));
}
double complex
ctan(double complex z)
{
/* ctan(z) = -I * ctanh(I * z) */
z = ctanh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
/* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */
z = ctanh(CMPLX(cimag(z), creal(z)));
return (CMPLX(cimag(z), creal(z)));
}

17
lib/msun/src/s_ctanhf.c
View File

@ -51,18 +51,19 @@ ctanhf(float complex z)
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
return (cpackf(x, (y == 0 ? y : x * y)));
return (CMPLXF((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
SET_FLOAT_WORD(x, hx - 0x40000000);
return (cpackf(x,
return (CMPLXF(x,
copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));
}
if (!isfinite(y))
return (cpackf(y - y, y - y));
return (CMPLXF(y - y, y - y));
if (ix >= 0x41300000) { /* x >= 11 */
if (ix >= 0x41300000) { /* |x| >= 11 */
float exp_mx = expf(-fabsf(x));
return (cpackf(copysignf(1, x),
return (CMPLXF(copysignf(1, x),
4 * sinf(y) * cosf(y) * exp_mx * exp_mx));
}
@ -71,14 +72,14 @@ ctanhf(float complex z)
s = sinhf(x);
rho = sqrtf(1 + s * s);
denom = 1 + beta * s * s;
return (cpackf((beta * rho * s) / denom, t / denom));
return (CMPLXF((beta * rho * s) / denom, t / denom));
}
float complex
ctanf(float complex z)
{
z = ctanhf(cpackf(-cimagf(z), crealf(z)));
return (cpackf(cimagf(z), -crealf(z)));
z = ctanhf(CMPLXF(cimagf(z), crealf(z)));
return (CMPLXF(cimagf(z), crealf(z)));
}

26
lib/msun/src/s_scalbln.c
View File

@ -27,32 +27,28 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <limits.h>
#include <math.h>
double
scalbln (double x, long n)
{
int in;
#define NMAX 65536
#define NMIN -65536
in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n;
return (scalbn(x, in));
double
scalbln(double x, long n)
{
return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}
float
scalblnf (float x, long n)
scalblnf(float x, long n)
{
int in;
in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n;
return (scalbnf(x, in));
return (scalbnf(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}
long double
scalblnl (long double x, long n)
scalblnl(long double x, long n)
{
int in;
in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n;
return (scalbnl(x, in));
return (scalbnl(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}