Improve cosmo's conformance to libc-test

This change addresses various open source compatibility issues, so that
we pass 313/411 of the tests in https://github.com/jart/libc-test where
earlier today we were passing about 30/411 of them, due to header toil.
Please note that Glibc only passes 341/411 so 313 today is pretty good!

- Make the conformance of libc/isystem/ headers nearly perfect
- Import more of the remaining math library routines from Musl
- Fix inconsistencies with type signatures of calls like umask
- Write tests for getpriority/setpriority which work great now
- conform to `struct sockaddr *` on remaining socket functions
- Import a bunch of uninteresting stdlib functions e.g. rand48
- Introduce readdir_r, scandir, pthread_kill, sigsetjmp, etc..

Follow the instructions in our `tool/scripts/cosmocc` toolchain to run
these tests yourself. You use `make CC=cosmocc` on the test repository
This commit is contained in:
Justine Tunney 2022-10-10 17:52:41 -07:00
parent 467a332e38
commit e557058ac8
No known key found for this signature in database
GPG key ID: BE714B4575D6E328
189 changed files with 5091 additions and 884 deletions

View file

@ -35,7 +35,6 @@ Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
/* clang-format off */
/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>

View file

@ -33,7 +33,7 @@ asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
/* clang-format off */
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */
/*

View file

@ -1,29 +1,148 @@
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi
Copyright 2021 Justine Alexandra Roberts Tunney
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
Musl Libc
Copyright © 2005-2020 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
fdlibm (fdlibm license)\\n\
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c */
/*
* ====================================================
* 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.
* ====================================================
*
*/
/* atan2(y,x)
* Method :
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
* 2. Reduce x to positive by (if x and y are unexceptional):
* ARG (x+iy) = arctan(y/x) ... if x > 0,
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
*
* Special cases:
*
* ATAN2((anything), NaN ) is NaN;
* ATAN2(NAN , (anything) ) is NaN;
* ATAN2(+-0, +(anything but NaN)) is +-0 ;
* ATAN2(+-0, -(anything but NaN)) is +-pi ;
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
* ATAN2(+-INF,+INF ) is +-pi/4 ;
* ATAN2(+-INF,-INF ) is +-3pi/4;
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const double
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
/**
* Returns arc tangent of 𝑦/𝑥.
* @note the greatest of all libm functions
*/
double atan2(double y, double x) {
long double st;
asm("fpatan" : "=t"(st) : "0"((long double)x), "u"((long double)y) : "st(1)");
return st;
double atan2(double y, double x)
{
double z;
uint32_t m,lx,ly,ix,iy;
if (isnan(x) || isnan(y))
return x+y;
EXTRACT_WORDS(ix, lx, x);
EXTRACT_WORDS(iy, ly, y);
if ((ix-0x3ff00000 | lx) == 0) /* x = 1.0 */
return atan(y);
m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */
ix = ix & 0x7fffffff;
iy = iy & 0x7fffffff;
/* when y = 0 */
if ((iy|ly) == 0) {
switch(m) {
case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */
case 2: return pi; /* atan(+0,-anything) = pi */
case 3: return -pi; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
if ((ix|lx) == 0)
return m&1 ? -pi/2 : pi/2;
/* when x is INF */
if (ix == 0x7ff00000) {
if (iy == 0x7ff00000) {
switch(m) {
case 0: return pi/4; /* atan(+INF,+INF) */
case 1: return -pi/4; /* atan(-INF,+INF) */
case 2: return 3*pi/4; /* atan(+INF,-INF) */
case 3: return -3*pi/4; /* atan(-INF,-INF) */
}
} else {
switch(m) {
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi; /* atan(+...,-INF) */
case 3: return -pi; /* atan(-...,-INF) */
}
}
}
/* |y/x| > 0x1p64 */
if (ix+(64<<20) < iy || iy == 0x7ff00000)
return m&1 ? -pi/2 : pi/2;
/* z = atan(|y/x|) without spurious underflow */
if ((m&2) && iy+(64<<20) < ix) /* |y/x| < 0x1p-64, x<0 */
z = 0;
else
z = atan(fabs(y/x));
switch (m) {
case 0: return z; /* atan(+,+) */
case 1: return -z; /* atan(-,+) */
case 2: return pi - (z-pi_lo); /* atan(+,-) */
default: /* case 3 */
return (z-pi_lo) - pi; /* atan(-,-) */
}
}

View file

@ -33,7 +33,7 @@ asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
/* clang-format off */
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */
/*

View file

@ -1,27 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
// Returns 𝑒^x-1.
//
// @param 𝑥 is double scalar in low half of %xmm0
// @return double scalar in low half of %xmm0
expm1: ezlea expm1l,ax
jmp _d2ld2
.endfn expm1,globl

238
libc/tinymath/expm1.c Normal file
View file

@ -0,0 +1,238 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/internal.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* expm1(x)
* Returns exp(x)-1, the exponential of x minus 1.
*
* Method
* 1. Argument reduction:
* Given x, find r and integer k such that
*
* x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
*
* Here a correction term c will be computed to compensate
* the error in r when rounded to a floating-point number.
*
* 2. Approximating expm1(r) by a special rational function on
* the interval [0,0.34658]:
* Since
* r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
* we define R1(r*r) by
* r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
* That is,
* R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
* = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
* = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
* We use a special Remez algorithm on [0,0.347] to generate
* a polynomial of degree 5 in r*r to approximate R1. The
* maximum error of this polynomial approximation is bounded
* by 2**-61. In other words,
* R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
* where Q1 = -1.6666666666666567384E-2,
* Q2 = 3.9682539681370365873E-4,
* Q3 = -9.9206344733435987357E-6,
* Q4 = 2.5051361420808517002E-7,
* Q5 = -6.2843505682382617102E-9;
* z = r*r,
* with error bounded by
* | 5 | -61
* | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
* | |
*
* expm1(r) = exp(r)-1 is then computed by the following
* specific way which minimize the accumulation rounding error:
* 2 3
* r r [ 3 - (R1 + R1*r/2) ]
* expm1(r) = r + --- + --- * [--------------------]
* 2 2 [ 6 - r*(3 - R1*r/2) ]
*
* To compensate the error in the argument reduction, we use
* expm1(r+c) = expm1(r) + c + expm1(r)*c
* ~ expm1(r) + c + r*c
* Thus c+r*c will be added in as the correction terms for
* expm1(r+c). Now rearrange the term to avoid optimization
* screw up:
* ( 2 2 )
* ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
* expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
* ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
* ( )
*
* = r - E
* 3. Scale back to obtain expm1(x):
* From step 1, we have
* expm1(x) = either 2^k*[expm1(r)+1] - 1
* = or 2^k*[expm1(r) + (1-2^-k)]
* 4. Implementation notes:
* (A). To save one multiplication, we scale the coefficient Qi
* to Qi*2^i, and replace z by (x^2)/2.
* (B). To achieve maximum accuracy, we compute expm1(x) by
* (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
* (ii) if k=0, return r-E
* (iii) if k=-1, return 0.5*(r-E)-0.5
* (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
* else return 1.0+2.0*(r-E);
* (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
* (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
* (vii) return 2^k(1-((E+2^-k)-r))
*
* Special cases:
* expm1(INF) is INF, expm1(NaN) is NaN;
* expm1(-INF) is -1, and
* for finite argument, only expm1(0)=0 is exact.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Misc. info.
* For IEEE double
* if x > 7.09782712893383973096e+02 then expm1(x) overflow
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const double
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
/**
* Returns 𝑒^𝑥-𝟷.
*/
double expm1(double x)
{
double_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
union {double f; uint64_t i;} u = {x};
uint32_t hx = u.i>>32 & 0x7fffffff;
int k, sign = u.i>>63;
/* filter out huge and non-finite argument */
if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */
if (isnan(x))
return x;
if (sign)
return -1;
if (x > o_threshold) {
x *= 0x1p1023;
return x;
}
}
/* argument reduction */
if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
if (!sign) {
hi = x - ln2_hi;
lo = ln2_lo;
k = 1;
} else {
hi = x + ln2_hi;
lo = -ln2_lo;
k = -1;
}
} else {
k = invln2*x + (sign ? -0.5 : 0.5);
t = k;
hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
lo = t*ln2_lo;
}
x = hi-lo;
c = (hi-x)-lo;
} else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */
if (hx < 0x00100000)
FORCE_EVAL((float)x);
return x;
} else
k = 0;
/* x is now in primary range */
hfx = 0.5*x;
hxs = x*hfx;
r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */
return x - (x*e-hxs);
e = x*(e-c) - c;
e -= hxs;
/* exp(x) ~ 2^k (x_reduced - e + 1) */
if (k == -1)
return 0.5*(x-e) - 0.5;
if (k == 1) {
if (x < -0.25)
return -2.0*(e-(x+0.5));
return 1.0+2.0*(x-e);
}
u.i = (uint64_t)(0x3ff + k)<<52; /* 2^k */
twopk = u.f;
if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */
y = x - e + 1.0;
if (k == 1024)
y = y*2.0*0x1p1023;
else
y = y*twopk;
return y - 1.0;
}
u.i = (uint64_t)(0x3ff - k)<<52; /* 2^-k */
if (k < 20)
y = (x-e+(1-u.f))*twopk;
else
y = (x-(e+u.f)+1)*twopk;
return y;
}

View file

@ -1,58 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
// Returns largest integral not greater than 𝑥.
//
// @param 𝑥 is double scalar in low half of %xmm0
// @return double scalar in low half of %xmm0
floor: .leafprologue
.profilable
movsd 4f(%rip),%xmm3
movsd 2f(%rip),%xmm4
movapd %xmm0,%xmm2
movapd %xmm0,%xmm1
andpd %xmm3,%xmm2
ucomisd %xmm2,%xmm4
jbe 1f
cvttsd2siq %xmm0,%rax
pxor %xmm2,%xmm2
movsd 3f(%rip),%xmm4
andnpd %xmm1,%xmm3
cvtsi2sdq %rax,%xmm2
movapd %xmm2,%xmm5
cmpnlesd %xmm0,%xmm5
movapd %xmm5,%xmm0
andpd %xmm4,%xmm0
subsd %xmm0,%xmm2
movapd %xmm2,%xmm0
orpd %xmm3,%xmm0
1: .leafepilogue
.endfn floor,globl
.rodata.cst8
2: .long 0x00000000
.long 0x43300000
3: .long 0x00000000
.long 0x3ff00000
.rodata.cst16
4: .long 0xffffffff
.long 0x7fffffff
.long 0x00000000
.long 0x00000000

65
libc/tinymath/floor.c Normal file
View file

@ -0,0 +1,65 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/internal.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
#define EPS DBL_EPSILON
#elif FLT_EVAL_METHOD==2
#define EPS LDBL_EPSILON
#endif
static const double_t toint = 1/EPS;
double floor(double x)
{
union {double f; uint64_t i;} u = {x};
int e = u.i >> 52 & 0x7ff;
double_t y;
if (e >= 0x3ff+52 || x == 0)
return x;
/* y = int(x) - x, where int(x) is an integer neighbor of x */
if (u.i >> 63)
y = x - toint + toint - x;
else
y = x + toint - toint - x;
/* special case because of non-nearest rounding modes */
if (e <= 0x3ff-1) {
FORCE_EVAL(y);
return u.i >> 63 ? -1 : 0;
}
if (y > 0)
return x + y - 1;
return x + y;
}

View file

@ -1,56 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
// Returns largest integral not greater than 𝑥.
//
// @param 𝑥 is float scalar in low quarter of %xmm0
// @return float scalar in low quarter of %xmm0
floorf: .leafprologue
.profilable
movss 4f(%rip),%xmm3
movss 2f(%rip),%xmm4
movaps %xmm0,%xmm2
movaps %xmm0,%xmm1
andps %xmm3,%xmm2
ucomiss %xmm2,%xmm4
jbe 1f
cvttss2sil %xmm0,%eax
pxor %xmm2,%xmm2
movss 3f(%rip),%xmm4
andnps %xmm1,%xmm3
cvtsi2ssl %eax,%xmm2
movaps %xmm2,%xmm5
cmpnless %xmm0,%xmm5
movaps %xmm5,%xmm0
andps %xmm4,%xmm0
subss %xmm0,%xmm2
movaps %xmm2,%xmm0
orps %xmm3,%xmm0
1: .leafepilogue
.endfn floorf,globl
.rodata.cst4
2: .long 0x4b000000
3: .long 0x3f800000
.rodata.cst16
4: .long 0x7fffffff
.long 0x00000000
.long 0x00000000
.long 0x00000000

61
libc/tinymath/floorf.c Normal file
View file

@ -0,0 +1,61 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/tinymath/internal.h"
#include "third_party/libcxx/math.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
float floorf(float x)
{
union {float f; uint32_t i;} u = {x};
int e = (int)(u.i >> 23 & 0xff) - 0x7f;
uint32_t m;
if (e >= 23)
return x;
if (e >= 0) {
m = 0x007fffff >> e;
if ((u.i & m) == 0)
return x;
FORCE_EVAL(x + 0x1p120f);
if (u.i >> 31)
u.i += m;
u.i &= ~m;
} else {
FORCE_EVAL(x + 0x1p120f);
if (u.i >> 31 == 0)
u.i = 0;
else if (u.i << 1)
u.f = -1.0;
}
return u.f;
}

View file

@ -25,7 +25,7 @@
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "third_party/libcxx/math.h"
#include "libc/math.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\

View file

@ -1,33 +1,101 @@
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi
Copyright 2021 Justine Alexandra Roberts Tunney
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/math.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64
#define SPLIT (0x1p32 + 1)
#else
#define SPLIT (0x1p27 + 1)
#endif
static void sq(double_t *hi, double_t *lo, double x)
{
double_t xh, xl, xc;
xc = (double_t)x*SPLIT;
xh = x - xc + xc;
xl = x - xh;
*hi = (double_t)x*x;
*lo = xh*xh - *hi + 2*xh*xl + xl*xl;
}
/**
* Returns euclidean distance.
*/
double hypot(double a, double b) {
double r, t;
if (isinf(a) || isinf(b)) return INFINITY;
a = fabs(a);
b = fabs(b);
if (a < b) t = b, b = a, a = t;
if (!a) return b;
r = b / a;
return a * sqrt(1 + r * r);
double hypot(double x, double y)
{
union {double f; uint64_t i;} ux = {x}, uy = {y}, ut;
int ex, ey;
double_t hx, lx, hy, ly, z;
/* arrange |x| >= |y| */
ux.i &= -1ULL>>1;
uy.i &= -1ULL>>1;
if (ux.i < uy.i) {
ut = ux;
ux = uy;
uy = ut;
}
/* special cases */
ex = ux.i>>52;
ey = uy.i>>52;
x = ux.f;
y = uy.f;
/* note: hypot(inf,nan) == inf */
if (ey == 0x7ff)
return y;
if (ex == 0x7ff || uy.i == 0)
return x;
/* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
/* 64 difference is enough for ld80 double_t */
if (ex - ey > 64)
return x + y;
/* precise sqrt argument in nearest rounding mode without overflow */
/* xh*xh must not overflow and xl*xl must not underflow in sq */
z = 1;
if (ex > 0x3ff+510) {
z = 0x1p700;
x *= 0x1p-700;
y *= 0x1p-700;
} else if (ey < 0x3ff-450) {
z = 0x1p-700;
x *= 0x1p700;
y *= 0x1p700;
}
sq(&hx, &lx, x);
sq(&hy, &ly, y);
return z*sqrt(ly+lx+hy+hx);
}

View file

@ -1,33 +1,67 @@
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi
Copyright 2021 Justine Alexandra Roberts Tunney
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/math.h"
/**
* Returns euclidean distance.
*/
float hypotf(float a, float b) {
float r, t;
if (isinf(a) || isinf(b)) return INFINITY;
a = fabsf(a);
b = fabsf(b);
if (a < b) t = b, b = a, a = t;
if (!a) return b;
r = b / a;
return a * sqrtf(1 + r * r);
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
float hypotf(float x, float y)
{
union {float f; uint32_t i;} ux = {x}, uy = {y}, ut;
float_t z;
ux.i &= -1U>>1;
uy.i &= -1U>>1;
if (ux.i < uy.i) {
ut = ux;
ux = uy;
uy = ut;
}
x = ux.f;
y = uy.f;
if (uy.i == 0xff<<23)
return y;
if (ux.i >= 0xff<<23 || uy.i == 0 || ux.i - uy.i >= 25<<23)
return x + y;
z = 1;
if (ux.i >= (0x7f+60)<<23) {
z = 0x1p90f;
x *= 0x1p-90f;
y *= 0x1p-90f;
} else if (uy.i < (0x7f-60)<<23) {
z = 0x1p-90f;
x *= 0x1p90f;
y *= 0x1p90f;
}
return z*sqrtf((double)x*x + (double)y*y);
}

View file

@ -1,33 +1,100 @@
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi
Copyright 2021 Justine Alexandra Roberts Tunney
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/ldshape.internal.h"
/**
* Returns euclidean distance.
*/
long double hypotl(long double a, long double b) {
long double r, t;
if (isinf(a) || isinf(b)) return INFINITY;
a = fabsl(a);
b = fabsl(b);
if (a < b) t = b, b = a, a = t;
if (!a) return b;
r = b / a;
return a * sqrtl(1 + r * r);
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double hypotl(long double x, long double y)
{
return hypot(x, y);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#if LDBL_MANT_DIG == 64
#define SPLIT (0x1p32L+1)
#elif LDBL_MANT_DIG == 113
#define SPLIT (0x1p57L+1)
#endif
static void sq(long double *hi, long double *lo, long double x)
{
long double xh, xl, xc;
xc = x*SPLIT;
xh = x - xc + xc;
xl = x - xh;
*hi = x*x;
*lo = xh*xh - *hi + 2*xh*xl + xl*xl;
}
long double hypotl(long double x, long double y)
{
union ldshape ux = {x}, uy = {y};
int ex, ey;
long double hx, lx, hy, ly, z;
ux.i.se &= 0x7fff;
uy.i.se &= 0x7fff;
if (ux.i.se < uy.i.se) {
ex = uy.i.se;
ey = ux.i.se;
x = uy.f;
y = ux.f;
} else {
ex = ux.i.se;
ey = uy.i.se;
x = ux.f;
y = uy.f;
}
if (ex == 0x7fff && isinf(y))
return y;
if (ex == 0x7fff || y == 0)
return x;
if (ex - ey > LDBL_MANT_DIG)
return x + y;
z = 1;
if (ex > 0x3fff+8000) {
z = 0x1p10000L;
x *= 0x1p-10000L;
y *= 0x1p-10000L;
} else if (ey < 0x3fff-8000) {
z = 0x1p-10000L;
x *= 0x1p10000L;
y *= 0x1p10000L;
}
sq(&hx, &lx, x);
sq(&hy, &ly, y);
return z*sqrtl(ly+lx+hy+hx);
}
#endif

409
libc/tinymath/j0.c Normal file
View file

@ -0,0 +1,409 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_j0.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/* j0(x), y0(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j0(x):
* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
* 2. Reduce x to |x| since j0(x)=j0(-x), and
* for x in (0,2)
* j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
* (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
* for x in (2,inf)
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* as follow:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (cos(x) + sin(x))
* sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j0(nan)= nan
* j0(0) = 1
* j0(inf) = 0
*
* Method -- y0(x):
* 1. For x<2.
* Since
* y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
* therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
* We use the following function to approximate y0,
* y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
* where
* U(z) = u00 + u01*z + ... + u06*z^6
* V(z) = 1 + v01*z + ... + v04*z^4
* with absolute approximation error bounded by 2**-72.
* Note: For tiny x, U/V = u0 and j0(x)~1, hence
* y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
* 2. For x>=2.
* y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* by the method mentioned above.
* 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
*/
static double pzero(double), qzero(double);
static const double
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
/* common method when |x|>=2 */
static double common(uint32_t ix, double x, int y0)
{
double s,c,ss,cc,z;
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x-pi/4)-q0(x)*sin(x-pi/4))
* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x-pi/4)+q0(x)*cos(x-pi/4))
*
* sin(x-pi/4) = (sin(x) - cos(x))/sqrt(2)
* cos(x-pi/4) = (sin(x) + cos(x))/sqrt(2)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
*/
s = sin(x);
c = cos(x);
if (y0)
c = -c;
cc = s+c;
/* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */
if (ix < 0x7fe00000) {
ss = s-c;
z = -cos(2*x);
if (s*c < 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x48000000) {
if (y0)
ss = -ss;
cc = pzero(x)*cc-qzero(x)*ss;
}
}
return invsqrtpi*cc/sqrt(x);
}
/* R0/S0 on [0, 2.00] */
static const double
R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
R05 = -4.61832688532103189199e-09, /* 0xBE33D5E7, 0x73D63FCE */
S01 = 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
S02 = 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
double j0(double x)
{
double z,r,s;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
/* j0(+-inf)=0, j0(nan)=nan */
if (ix >= 0x7ff00000)
return 1/(x*x);
x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2 */
/* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */
return common(ix,x,0);
}
/* 1 - x*x/4 + x*x*R(x^2)/S(x^2) */
if (ix >= 0x3f200000) { /* |x| >= 2**-13 */
/* up to 4ulp error close to 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1+z*(S01+z*(S02+z*(S03+z*S04)));
return (1+x/2)*(1-x/2) + z*(r/s);
}
/* 1 - x*x/4 */
/* prevent underflow */
/* inexact should be raised when x!=0, this is not done correctly */
if (ix >= 0x38000000) /* |x| >= 2**-127 */
x = 0.25*x*x;
return 1 - x;
}
static const double
u00 = -7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
u01 = 1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
u02 = -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
u03 = 3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
u04 = -3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
u05 = 1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
u06 = -3.98205194132103398453e-11, /* 0xBDC5E43D, 0x693FB3C8 */
v01 = 1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
v02 = 7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
v03 = 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
double y0(double x)
{
double z,u,v;
uint32_t ix,lx;
EXTRACT_WORDS(ix, lx, x);
/* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */
if ((ix<<1 | lx) == 0)
return -1/0.0;
if (ix>>31)
return 0/0.0;
if (ix >= 0x7ff00000)
return 1/x;
if (ix >= 0x40000000) { /* x >= 2 */
/* large ulp errors near zeros: 3.958, 7.086,.. */
return common(ix,x,1);
}
/* U(x^2)/V(x^2) + (2/pi)*j0(x)*log(x) */
if (ix >= 0x3e400000) { /* x >= 2**-27 */
/* large ulp error near the first zero, x ~= 0.89 */
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x));
}
return u00 + tpi*log(x);
}
/* The asymptotic expansions of pzero is
* 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
* For x >= 2, We approximate pzero by
* pzero(x) = 1 + (R/S)
* where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
* S = 1 + pS0*s^2 + ... + pS4*s^10
* and
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
*/
static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
-8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
-2.57063105679704847262e+02, /* 0xC0701102, 0x7B19E863 */
-2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
-5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
};
static const double pS8[5] = {
1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
1.16752972564375915681e+05, /* 0x40FC810F, 0x8F9FA9BD */
4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
};
static const double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
-7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
-4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
-6.76747652265167261021e+01, /* 0xC050EB2F, 0x5A7D1783 */
-3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
-3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
};
static const double pS5[5] = {
6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
9.62544514357774460223e+03, /* 0x40C2CCB8, 0xFA76FA38 */
2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
};
static const double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
-7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
-2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
-2.19659774734883086467e+01, /* 0xC035F74A, 0x4CB94E14 */
-5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
-3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
};
static const double pS3[5] = {
3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
1.12799679856907414432e+03, /* 0x40919FFC, 0xB8C39B7E */
1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
};
static const double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
-7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
-1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
-7.63569613823527770791e+00, /* 0xC01E8AF3, 0xEDAFA7F3 */
-1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
-3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
};
static const double pS2[5] = {
2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
};
static double pzero(double x)
{
const double *p,*q;
double_t z,r,s;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = pR8; q = pS8;}
else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0 + r/s;
}
/* For x >= 8, the asymptotic expansions of qzero is
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
* We approximate pzero by
* qzero(x) = s*(-1.25 + (R/S))
* where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
* S = 1 + qS0*s^2 + ... + qS5*s^12
* and
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
*/
static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
5.57673380256401856059e+02, /* 0x40816D63, 0x15301825 */
8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
};
static const double qS8[6] = {
1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
8.03309257119514397345e+05, /* 0x412883DA, 0x83A52B43 */
8.40501579819060512818e+05, /* 0x4129A66B, 0x28DE0B3D */
-3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
};
static const double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
1.35111577286449829671e+02, /* 0x4060E392, 0x0A8788E9 */
1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
};
static const double qS5[6] = {
8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
5.67511122894947329769e+04, /* 0x40EBB5E3, 0x97E02372 */
3.59767538425114471465e+04, /* 0x40E19118, 0x1F7A54A0 */
-5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
};
static const double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
4.26218440745412650017e+01, /* 0x40454F98, 0x962DAEDD */
1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
};
static const double qS3[6] = {
4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
6.46042516752568917582e+03, /* 0x40B93C6C, 0xD7C76A28 */
2.51633368920368957333e+03, /* 0x40A3A8AA, 0xD94FB1C0 */
-1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
};
static const double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
1.44956029347885735348e+01, /* 0x402CFDBF, 0xAAF96FE5 */
3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
};
static const double qS2[6] = {
3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
-5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
};
static double qzero(double x)
{
const double *p,*q;
double_t s,r,z;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = qR8; q = qS8;}
else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125 + r/s)/x;
}

347
libc/tinymath/j0f.c Normal file
View file

@ -0,0 +1,347 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
static float pzerof(float), qzerof(float);
static const float
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01; /* 0x3f22f983 */
static float common(uint32_t ix, float x, int y0)
{
float z,s,c,ss,cc;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
s = sinf(x);
c = cosf(x);
if (y0)
c = -c;
cc = s+c;
if (ix < 0x7f000000) {
ss = s-c;
z = -cosf(2*x);
if (s*c < 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x58800000) {
if (y0)
ss = -ss;
cc = pzerof(x)*cc-qzerof(x)*ss;
}
}
return invsqrtpi*cc/sqrtf(x);
}
/* R0/S0 on [0, 2.00] */
static const float
R02 = 1.5625000000e-02, /* 0x3c800000 */
R03 = -1.8997929874e-04, /* 0xb947352e */
R04 = 1.8295404516e-06, /* 0x35f58e88 */
R05 = -4.6183270541e-09, /* 0xb19eaf3c */
S01 = 1.5619102865e-02, /* 0x3c7fe744 */
S02 = 1.1692678527e-04, /* 0x38f53697 */
S03 = 5.1354652442e-07, /* 0x3509daa6 */
S04 = 1.1661400734e-09; /* 0x30a045e8 */
float j0f(float x)
{
float z,r,s;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x7f800000)
return 1/(x*x);
x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2 */
/* large ulp error near zeros */
return common(ix, x, 0);
}
if (ix >= 0x3a000000) { /* |x| >= 2**-11 */
/* up to 4ulp error near 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1+z*(S01+z*(S02+z*(S03+z*S04)));
return (1+x/2)*(1-x/2) + z*(r/s);
}
if (ix >= 0x21800000) /* |x| >= 2**-60 */
x = 0.25f*x*x;
return 1 - x;
}
static const float
u00 = -7.3804296553e-02, /* 0xbd9726b5 */
u01 = 1.7666645348e-01, /* 0x3e34e80d */
u02 = -1.3818567619e-02, /* 0xbc626746 */
u03 = 3.4745343146e-04, /* 0x39b62a69 */
u04 = -3.8140706238e-06, /* 0xb67ff53c */
u05 = 1.9559013964e-08, /* 0x32a802ba */
u06 = -3.9820518410e-11, /* 0xae2f21eb */
v01 = 1.2730483897e-02, /* 0x3c509385 */
v02 = 7.6006865129e-05, /* 0x389f65e0 */
v03 = 2.5915085189e-07, /* 0x348b216c */
v04 = 4.4111031494e-10; /* 0x2ff280c2 */
float y0f(float x)
{
float z,u,v;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
if ((ix & 0x7fffffff) == 0)
return -1/0.0f;
if (ix>>31)
return 0/0.0f;
if (ix >= 0x7f800000)
return 1/x;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* large ulp error near zeros */
return common(ix,x,1);
}
if (ix >= 0x39000000) { /* x >= 2**-13 */
/* large ulp error at x ~= 0.89 */
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x));
}
return u00 + tpi*logf(x);
}
/* The asymptotic expansions of pzero is
* 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
* For x >= 2, We approximate pzero by
* pzero(x) = 1 + (R/S)
* where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
* S = 1 + pS0*s^2 + ... + pS4*s^10
* and
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
*/
static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-7.0312500000e-02, /* 0xbd900000 */
-8.0816707611e+00, /* 0xc1014e86 */
-2.5706311035e+02, /* 0xc3808814 */
-2.4852163086e+03, /* 0xc51b5376 */
-5.2530439453e+03, /* 0xc5a4285a */
};
static const float pS8[5] = {
1.1653436279e+02, /* 0x42e91198 */
3.8337448730e+03, /* 0x456f9beb */
4.0597855469e+04, /* 0x471e95db */
1.1675296875e+05, /* 0x47e4087c */
4.7627726562e+04, /* 0x473a0bba */
};
static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.1412546255e-11, /* 0xad48c58a */
-7.0312492549e-02, /* 0xbd8fffff */
-4.1596107483e+00, /* 0xc0851b88 */
-6.7674766541e+01, /* 0xc287597b */
-3.3123129272e+02, /* 0xc3a59d9b */
-3.4643338013e+02, /* 0xc3ad3779 */
};
static const float pS5[5] = {
6.0753936768e+01, /* 0x42730408 */
1.0512523193e+03, /* 0x44836813 */
5.9789707031e+03, /* 0x45bad7c4 */
9.6254453125e+03, /* 0x461665c8 */
2.4060581055e+03, /* 0x451660ee */
};
static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.5470459075e-09, /* 0xb12f081b */
-7.0311963558e-02, /* 0xbd8fffb8 */
-2.4090321064e+00, /* 0xc01a2d95 */
-2.1965976715e+01, /* 0xc1afba52 */
-5.8079170227e+01, /* 0xc2685112 */
-3.1447946548e+01, /* 0xc1fb9565 */
};
static const float pS3[5] = {
3.5856033325e+01, /* 0x420f6c94 */
3.6151397705e+02, /* 0x43b4c1ca */
1.1936077881e+03, /* 0x44953373 */
1.1279968262e+03, /* 0x448cffe6 */
1.7358093262e+02, /* 0x432d94b8 */
};
static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.8753431271e-08, /* 0xb3be98b7 */
-7.0303097367e-02, /* 0xbd8ffb12 */
-1.4507384300e+00, /* 0xbfb9b1cc */
-7.6356959343e+00, /* 0xc0f4579f */
-1.1193166733e+01, /* 0xc1331736 */
-3.2336456776e+00, /* 0xc04ef40d */
};
static const float pS2[5] = {
2.2220300674e+01, /* 0x41b1c32d */
1.3620678711e+02, /* 0x430834f0 */
2.7047027588e+02, /* 0x43873c32 */
1.5387539673e+02, /* 0x4319e01a */
1.4657617569e+01, /* 0x416a859a */
};
static float pzerof(float x)
{
const float *p,*q;
float_t z,r,s;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = pR8; q = pS8;}
else if (ix >= 0x409173eb){p = pR5; q = pS5;}
else if (ix >= 0x4036d917){p = pR3; q = pS3;}
else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0f + r/s;
}
/* For x >= 8, the asymptotic expansions of qzero is
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
* We approximate pzero by
* qzero(x) = s*(-1.25 + (R/S))
* where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
* S = 1 + qS0*s^2 + ... + qS5*s^12
* and
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
*/
static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
7.3242187500e-02, /* 0x3d960000 */
1.1768206596e+01, /* 0x413c4a93 */
5.5767340088e+02, /* 0x440b6b19 */
8.8591972656e+03, /* 0x460a6cca */
3.7014625000e+04, /* 0x471096a0 */
};
static const float qS8[6] = {
1.6377603149e+02, /* 0x4323c6aa */
8.0983447266e+03, /* 0x45fd12c2 */
1.4253829688e+05, /* 0x480b3293 */
8.0330925000e+05, /* 0x49441ed4 */
8.4050156250e+05, /* 0x494d3359 */
-3.4389928125e+05, /* 0xc8a7eb69 */
};
static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.8408595828e-11, /* 0x2da1ec79 */
7.3242180049e-02, /* 0x3d95ffff */
5.8356351852e+00, /* 0x40babd86 */
1.3511157227e+02, /* 0x43071c90 */
1.0272437744e+03, /* 0x448067cd */
1.9899779053e+03, /* 0x44f8bf4b */
};
static const float qS5[6] = {
8.2776611328e+01, /* 0x42a58da0 */
2.0778142090e+03, /* 0x4501dd07 */
1.8847289062e+04, /* 0x46933e94 */
5.6751113281e+04, /* 0x475daf1d */
3.5976753906e+04, /* 0x470c88c1 */
-5.3543427734e+03, /* 0xc5a752be */
};
static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.3774099900e-09, /* 0x3196681b */
7.3241114616e-02, /* 0x3d95ff70 */
3.3442313671e+00, /* 0x405607e3 */
4.2621845245e+01, /* 0x422a7cc5 */
1.7080809021e+02, /* 0x432acedf */
1.6673394775e+02, /* 0x4326bbe4 */
};
static const float qS3[6] = {
4.8758872986e+01, /* 0x42430916 */
7.0968920898e+02, /* 0x44316c1c */
3.7041481934e+03, /* 0x4567825f */
6.4604252930e+03, /* 0x45c9e367 */
2.5163337402e+03, /* 0x451d4557 */
-1.4924745178e+02, /* 0xc3153f59 */
};
static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.5044444979e-07, /* 0x342189db */
7.3223426938e-02, /* 0x3d95f62a */
1.9981917143e+00, /* 0x3fffc4bf */
1.4495602608e+01, /* 0x4167edfd */
3.1666231155e+01, /* 0x41fd5471 */
1.6252708435e+01, /* 0x4182058c */
};
static const float qS2[6] = {
3.0365585327e+01, /* 0x41f2ecb8 */
2.6934811401e+02, /* 0x4386ac8f */
8.4478375244e+02, /* 0x44533229 */
8.8293585205e+02, /* 0x445cbbe5 */
2.1266638184e+02, /* 0x4354aa98 */
-5.3109550476e+00, /* 0xc0a9f358 */
};
static float qzerof(float x)
{
const float *p,*q;
float_t s,r,z;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = qR8; q = qS8;}
else if (ix >= 0x409173eb){p = qR5; q = qS5;}
else if (ix >= 0x4036d917){p = qR3; q = qS3;}
else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125f + r/s)/x;
}

396
libc/tinymath/j1.c Normal file
View file

@ -0,0 +1,396 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_j1.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/* j1(x), y1(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j1(x):
* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
* 2. Reduce x to |x| since j1(x)=-j1(-x), and
* for x in (0,2)
* j1(x) = x/2 + x*z*R0/S0, where z = x*x;
* (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 )
* for x in (2,inf)
* j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* as follow:
* cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (sin(x) + cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j1(nan)= nan
* j1(0) = 0
* j1(inf) = 0
*
* Method -- y1(x):
* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
* 2. For x<2.
* Since
* y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
* therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
* We use the following function to approximate y1,
* y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
* where for x in [0,2] (abs err less than 2**-65.89)
* U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
* V(z) = 1 + v0[0]*z + ... + v0[4]*z^5
* Note: For tiny x, 1/x dominate y1 and hence
* y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
* 3. For x>=2.
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* by method mentioned above.
*/
static double pone(double), qone(double);
static const double
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
static double common(uint32_t ix, double x, int y1, int sign)
{
double z,s,c,ss,cc;
/*
* j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x-3pi/4)-q1(x)*sin(x-3pi/4))
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x-3pi/4)+q1(x)*cos(x-3pi/4))
*
* sin(x-3pi/4) = -(sin(x) + cos(x))/sqrt(2)
* cos(x-3pi/4) = (sin(x) - cos(x))/sqrt(2)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
*/
s = sin(x);
if (y1)
s = -s;
c = cos(x);
cc = s-c;
if (ix < 0x7fe00000) {
/* avoid overflow in 2*x */
ss = -s-c;
z = cos(2*x);
if (s*c > 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x48000000) {
if (y1)
ss = -ss;
cc = pone(x)*cc-qone(x)*ss;
}
}
if (sign)
cc = -cc;
return invsqrtpi*cc/sqrt(x);
}
/* R0/S0 on [0,2] */
static const double
r00 = -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
r01 = 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
r02 = -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
r03 = 4.96727999609584448412e-08, /* 0x3E6AAAFA, 0x46CA0BD9 */
s01 = 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
s02 = 1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
s03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
double j1(double x)
{
double z,r,s;
uint32_t ix;
int sign;
GET_HIGH_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7ff00000)
return 1/(x*x);
if (ix >= 0x40000000) /* |x| >= 2 */
return common(ix, fabs(x), 0, sign);
if (ix >= 0x38000000) { /* |x| >= 2**-127 */
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
z = r/s;
} else
/* avoid underflow, raise inexact if x!=0 */
z = x;
return (0.5 + z)*x;
}
static const double U0[5] = {
-1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
-1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
-9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
};
static const double V0[5] = {
1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
};
double y1(double x)
{
double z,u,v;
uint32_t ix,lx;
EXTRACT_WORDS(ix, lx, x);
/* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */
if ((ix<<1 | lx) == 0)
return -1/0.0;
if (ix>>31)
return 0/0.0;
if (ix >= 0x7ff00000)
return 1/x;
if (ix >= 0x40000000) /* x >= 2 */
return common(ix, x, 1, 0);
if (ix < 0x3c900000) /* x < 2**-54 */
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = 1+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1(x)*log(x)-1/x);
}
/* For x >= 8, the asymptotic expansions of pone is
* 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
* We approximate pone by
* pone(x) = 1 + (R/S)
* where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
* S = 1 + ps0*s^2 + ... + ps4*s^10
* and
* | pone(x)-1-R/S | <= 2 ** ( -60.06)
*/
static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
};
static const double ps8[5] = {
1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
};
static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
};
static const double ps5[5] = {
5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
};
static const double pr3[6] = {
3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
};
static const double ps3[5] = {
3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
};
static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
};
static const double ps2[5] = {
2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
};
static double pone(double x)
{
const double *p,*q;
double_t z,r,s;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = pr8; q = ps8;}
else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
else /*ix >= 0x40000000*/ {p = pr2; q = ps2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0+ r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
* We approximate pone by
* qone(x) = s*(0.375 + (R/S))
* where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
* S = 1 + qs1*s^2 + ... + qs6*s^12
* and
* | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13)
*/
static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
-1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
-7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
-1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
-4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
};
static const double qs8[6] = {
1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
-2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
};
static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
-1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
-8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
-1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
-1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
-2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
};
static const double qs5[6] = {
8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
-4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
};
static const double qr3[6] = {
-5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
-1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
-4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
-5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
-2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
-2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
};
static const double qs3[6] = {
4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
-1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
};
static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
-1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
-2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
-1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
-4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
-2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
};
static const double qs2[6] = {
2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
-4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
};
static double qone(double x)
{
const double *p,*q;
double_t s,r,z;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = qr8; q = qs8;}
else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
else /*ix >= 0x40000000*/ {p = qr2; q = qs2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375 + r/s)/x;
}

343
libc/tinymath/j1f.c Normal file
View file

@ -0,0 +1,343 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_j1f.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
static float ponef(float), qonef(float);
static const float
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01; /* 0x3f22f983 */
static float common(uint32_t ix, float x, int y1, int sign)
{
double z,s,c,ss,cc;
s = sinf(x);
if (y1)
s = -s;
c = cosf(x);
cc = s-c;
if (ix < 0x7f000000) {
ss = -s-c;
z = cosf(2*x);
if (s*c > 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x58800000) {
if (y1)
ss = -ss;
cc = ponef(x)*cc-qonef(x)*ss;
}
}
if (sign)
cc = -cc;
return invsqrtpi*cc/sqrtf(x);
}
/* R0/S0 on [0,2] */
static const float
r00 = -6.2500000000e-02, /* 0xbd800000 */
r01 = 1.4070566976e-03, /* 0x3ab86cfd */
r02 = -1.5995563444e-05, /* 0xb7862e36 */
r03 = 4.9672799207e-08, /* 0x335557d2 */
s01 = 1.9153760746e-02, /* 0x3c9ce859 */
s02 = 1.8594678841e-04, /* 0x3942fab6 */
s03 = 1.1771846857e-06, /* 0x359dffc2 */
s04 = 5.0463624390e-09, /* 0x31ad6446 */
s05 = 1.2354227016e-11; /* 0x2d59567e */
float j1f(float x)
{
float z,r,s;
uint32_t ix;
int sign;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x7f800000)
return 1/(x*x);
if (ix >= 0x40000000) /* |x| >= 2 */
return common(ix, fabsf(x), 0, sign);
if (ix >= 0x39000000) { /* |x| >= 2**-13 */
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
z = 0.5f + r/s;
} else
z = 0.5f;
return z*x;
}
static const float U0[5] = {
-1.9605709612e-01, /* 0xbe48c331 */
5.0443872809e-02, /* 0x3d4e9e3c */
-1.9125689287e-03, /* 0xbafaaf2a */
2.3525259166e-05, /* 0x37c5581c */
-9.1909917899e-08, /* 0xb3c56003 */
};
static const float V0[5] = {
1.9916731864e-02, /* 0x3ca3286a */
2.0255257550e-04, /* 0x3954644b */
1.3560879779e-06, /* 0x35b602d4 */
6.2274145840e-09, /* 0x31d5f8eb */
1.6655924903e-11, /* 0x2d9281cf */
};
float y1f(float x)
{
float z,u,v;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
if ((ix & 0x7fffffff) == 0)
return -1/0.0f;
if (ix>>31)
return 0/0.0f;
if (ix >= 0x7f800000)
return 1/x;
if (ix >= 0x40000000) /* |x| >= 2.0 */
return common(ix,x,1,0);
if (ix < 0x33000000) /* x < 2**-25 */
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x);
}
/* For x >= 8, the asymptotic expansions of pone is
* 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
* We approximate pone by
* pone(x) = 1 + (R/S)
* where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
* S = 1 + ps0*s^2 + ... + ps4*s^10
* and
* | pone(x)-1-R/S | <= 2 ** ( -60.06)
*/
static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
1.1718750000e-01, /* 0x3df00000 */
1.3239480972e+01, /* 0x4153d4ea */
4.1205184937e+02, /* 0x43ce06a3 */
3.8747453613e+03, /* 0x45722bed */
7.9144794922e+03, /* 0x45f753d6 */
};
static const float ps8[5] = {
1.1420736694e+02, /* 0x42e46a2c */
3.6509309082e+03, /* 0x45642ee5 */
3.6956207031e+04, /* 0x47105c35 */
9.7602796875e+04, /* 0x47bea166 */
3.0804271484e+04, /* 0x46f0a88b */
};
static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.3199052094e-11, /* 0x2d68333f */
1.1718749255e-01, /* 0x3defffff */
6.8027510643e+00, /* 0x40d9b023 */
1.0830818176e+02, /* 0x42d89dca */
5.1763616943e+02, /* 0x440168b7 */
5.2871520996e+02, /* 0x44042dc6 */
};
static const float ps5[5] = {
5.9280597687e+01, /* 0x426d1f55 */
9.9140142822e+02, /* 0x4477d9b1 */
5.3532670898e+03, /* 0x45a74a23 */
7.8446904297e+03, /* 0x45f52586 */
1.5040468750e+03, /* 0x44bc0180 */
};
static const float pr3[6] = {
3.0250391081e-09, /* 0x314fe10d */
1.1718686670e-01, /* 0x3defffab */
3.9329774380e+00, /* 0x407bb5e7 */
3.5119403839e+01, /* 0x420c7a45 */
9.1055007935e+01, /* 0x42b61c2a */
4.8559066772e+01, /* 0x42423c7c */
};
static const float ps3[5] = {
3.4791309357e+01, /* 0x420b2a4d */
3.3676245117e+02, /* 0x43a86198 */
1.0468714600e+03, /* 0x4482dbe3 */
8.9081134033e+02, /* 0x445eb3ed */
1.0378793335e+02, /* 0x42cf936c */
};
static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.0771083225e-07, /* 0x33e74ea8 */
1.1717621982e-01, /* 0x3deffa16 */
2.3685150146e+00, /* 0x401795c0 */
1.2242610931e+01, /* 0x4143e1bc */
1.7693971634e+01, /* 0x418d8d41 */
5.0735230446e+00, /* 0x40a25a4d */
};
static const float ps2[5] = {
2.1436485291e+01, /* 0x41ab7dec */
1.2529022980e+02, /* 0x42fa9499 */
2.3227647400e+02, /* 0x436846c7 */
1.1767937469e+02, /* 0x42eb5bd7 */
8.3646392822e+00, /* 0x4105d590 */
};
static float ponef(float x)
{
const float *p,*q;
float_t z,r,s;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = pr8; q = ps8;}
else if (ix >= 0x409173eb){p = pr5; q = ps5;}
else if (ix >= 0x4036d917){p = pr3; q = ps3;}
else /*ix >= 0x40000000*/ {p = pr2; q = ps2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0f + r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
* We approximate pone by
* qone(x) = s*(0.375 + (R/S))
* where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
* S = 1 + qs1*s^2 + ... + qs6*s^12
* and
* | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13)
*/
static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-1.0253906250e-01, /* 0xbdd20000 */
-1.6271753311e+01, /* 0xc1822c8d */
-7.5960174561e+02, /* 0xc43de683 */
-1.1849806641e+04, /* 0xc639273a */
-4.8438511719e+04, /* 0xc73d3683 */
};
static const float qs8[6] = {
1.6139537048e+02, /* 0x43216537 */
7.8253862305e+03, /* 0x45f48b17 */
1.3387534375e+05, /* 0x4802bcd6 */
7.1965775000e+05, /* 0x492fb29c */
6.6660125000e+05, /* 0x4922be94 */
-2.9449025000e+05, /* 0xc88fcb48 */
};
static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.0897993405e-11, /* 0xadb7d219 */
-1.0253904760e-01, /* 0xbdd1fffe */
-8.0564479828e+00, /* 0xc100e736 */
-1.8366960144e+02, /* 0xc337ab6b */
-1.3731937256e+03, /* 0xc4aba633 */
-2.6124443359e+03, /* 0xc523471c */
};
static const float qs5[6] = {
8.1276550293e+01, /* 0x42a28d98 */
1.9917987061e+03, /* 0x44f8f98f */
1.7468484375e+04, /* 0x468878f8 */
4.9851425781e+04, /* 0x4742bb6d */
2.7948074219e+04, /* 0x46da5826 */
-4.7191835938e+03, /* 0xc5937978 */
};
static const float qr3[6] = {
-5.0783124372e-09, /* 0xb1ae7d4f */
-1.0253783315e-01, /* 0xbdd1ff5b */
-4.6101160049e+00, /* 0xc0938612 */
-5.7847221375e+01, /* 0xc267638e */
-2.2824453735e+02, /* 0xc3643e9a */
-2.1921012878e+02, /* 0xc35b35cb */
};
static const float qs3[6] = {
4.7665153503e+01, /* 0x423ea91e */
6.7386511230e+02, /* 0x4428775e */
3.3801528320e+03, /* 0x45534272 */
5.5477290039e+03, /* 0x45ad5dd5 */
1.9031191406e+03, /* 0x44ede3d0 */
-1.3520118713e+02, /* 0xc3073381 */
};
static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.7838172539e-07, /* 0xb43f8932 */
-1.0251704603e-01, /* 0xbdd1f475 */
-2.7522056103e+00, /* 0xc0302423 */
-1.9663616180e+01, /* 0xc19d4f16 */
-4.2325313568e+01, /* 0xc2294d1f */
-2.1371921539e+01, /* 0xc1aaf9b2 */
};
static const float qs2[6] = {
2.9533363342e+01, /* 0x41ec4454 */
2.5298155212e+02, /* 0x437cfb47 */
7.5750280762e+02, /* 0x443d602e */
7.3939318848e+02, /* 0x4438d92a */
1.5594900513e+02, /* 0x431bf2f2 */
-4.9594988823e+00, /* 0xc09eb437 */
};
static float qonef(float x)
{
const float *p,*q;
float_t s,r,z;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = qr8; q = qs8;}
else if (ix >= 0x409173eb){p = qr5; q = qs5;}
else if (ix >= 0x4036d917){p = qr3; q = qs3;}
else /*ix >= 0x40000000*/ {p = qr2; q = qs2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375f + r/s)/x;
}

314
libc/tinymath/jn.c Normal file
View file

@ -0,0 +1,314 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_jn.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* jn(n, x), yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
*
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
* Note 2. About jn(n,x), yn(n,x)
* For n=0, j0(x) is called,
* for n=1, j1(x) is called,
* for n<=x, forward recursion is used starting
* from values of j0(x) and j1(x).
* for n>x, a continued fraction approximation to
* j(n,x)/j(n-1,x) is evaluated and then backward
* recursion is used starting from a supposed value
* for j(n,x). The resulting value of j(0,x) is
* compared with the actual value to correct the
* supposed value of j(n,x).
*
* yn(n,x) is similar in all respects, except
* that forward recursion is used for all
* values of n>1.
*/
static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
double jn(int n, double x)
{
uint32_t ix, lx;
int nm1, i, sign;
double a, b, temp;
EXTRACT_WORDS(ix, lx, x);
sign = ix>>31;
ix &= 0x7fffffff;
if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */
return x;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
/* nm1 = |n|-1 is used instead of |n| to handle n==INT_MIN */
if (n == 0)
return j0(x);
if (n < 0) {
nm1 = -(n+1);
x = -x;
sign ^= 1;
} else
nm1 = n-1;
if (nm1 == 0)
return j1(x);
sign &= n; /* even n: 0, odd n: signbit(x) */
x = fabs(x);
if ((ix|lx) == 0 || ix == 0x7ff00000) /* if x is 0 or inf */
b = 0.0;
else if (nm1 < x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x52d00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch(nm1&3) {
case 0: temp = -cos(x)+sin(x); break;
case 1: temp = -cos(x)-sin(x); break;
case 2: temp = cos(x)-sin(x); break;
default:
case 3: temp = cos(x)+sin(x); break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = j0(x);
b = j1(x);
for (i=0; i<nm1; ) {
i++;
temp = b;
b = b*(2.0*i/x) - a; /* avoid underflow */
a = temp;
}
}
} else {
if (ix < 0x3e100000) { /* x < 2**-29 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (nm1 > 32) /* underflow */
b = 0.0;
else {
temp = x*0.5;
b = temp;
a = 1.0;
for (i=2; i<=nm1+1; i++) {
a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b = b/a;
}
} else {
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
/* determine k */
double t,q0,q1,w,h,z,tmp,nf;
int k;
nf = nm1 + 1.0;
w = 2*nf/x;
h = 2/x;
z = w+h;
q0 = w;
q1 = w*z - 1.0;
k = 1;
while (q1 < 1.0e9) {
k += 1;
z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
for (t=0.0, i=k; i>=0; i--)
t = 1/(2*(i+nf)/x - t);
a = t;
b = 1.0;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = nf*log(fabs(w));
if (tmp < 7.09782712893383973096e+02) {
for (i=nm1; i>0; i--) {
temp = b;
b = b*(2.0*i)/x - a;
a = temp;
}
} else {
for (i=nm1; i>0; i--) {
temp = b;
b = b*(2.0*i)/x - a;
a = temp;
/* scale b to avoid spurious overflow */
if (b > 0x1p500) {
a /= b;
t /= b;
b = 1.0;
}
}
}
z = j0(x);
w = j1(x);
if (fabs(z) >= fabs(w))
b = t*z/b;
else
b = t*w/a;
}
}
return sign ? -b : b;
}
double yn(int n, double x)
{
uint32_t ix, lx, ib;
int nm1, sign, i;
double a, b, temp;
EXTRACT_WORDS(ix, lx, x);
sign = ix>>31;
ix &= 0x7fffffff;
if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */
return x;
if (sign && (ix|lx)!=0) /* x < 0 */
return 0/0.0;
if (ix == 0x7ff00000)
return 0.0;
if (n == 0)
return y0(x);
if (n < 0) {
nm1 = -(n+1);
sign = n&1;
} else {
nm1 = n-1;
sign = 0;
}
if (nm1 == 0)
return sign ? -y1(x) : y1(x);
if (ix >= 0x52d00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch(nm1&3) {
case 0: temp = -sin(x)-cos(x); break;
case 1: temp = -sin(x)+cos(x); break;
case 2: temp = sin(x)+cos(x); break;
default:
case 3: temp = sin(x)-cos(x); break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = y0(x);
b = y1(x);
/* quit if b is -inf */
GET_HIGH_WORD(ib, b);
for (i=0; i<nm1 && ib!=0xfff00000; ){
i++;
temp = b;
b = (2.0*i/x)*b - a;
GET_HIGH_WORD(ib, b);
a = temp;
}
}
return sign ? -b : b;
}

235
libc/tinymath/jnf.c Normal file
View file

@ -0,0 +1,235 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_jnf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
float jnf(int n, float x)
{
uint32_t ix;
int nm1, sign, i;
float a, b, temp;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix > 0x7f800000) /* nan */
return x;
/* J(-n,x) = J(n,-x), use |n|-1 to avoid overflow in -n */
if (n == 0)
return j0f(x);
if (n < 0) {
nm1 = -(n+1);
x = -x;
sign ^= 1;
} else
nm1 = n-1;
if (nm1 == 0)
return j1f(x);
sign &= n; /* even n: 0, odd n: signbit(x) */
x = fabsf(x);
if (ix == 0 || ix == 0x7f800000) /* if x is 0 or inf */
b = 0.0f;
else if (nm1 < x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
b = j1f(x);
for (i=0; i<nm1; ){
i++;
temp = b;
b = b*(2.0f*i/x) - a;
a = temp;
}
} else {
if (ix < 0x35800000) { /* x < 2**-20 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (nm1 > 8) /* underflow */
nm1 = 8;
temp = 0.5f * x;
b = temp;
a = 1.0f;
for (i=2; i<=nm1+1; i++) {
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b = b/a;
} else {
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
/* determine k */
float t,q0,q1,w,h,z,tmp,nf;
int k;
nf = nm1+1.0f;
w = 2*nf/x;
h = 2/x;
z = w+h;
q0 = w;
q1 = w*z - 1.0f;
k = 1;
while (q1 < 1.0e4f) {
k += 1;
z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
for (t=0.0f, i=k; i>=0; i--)
t = 1.0f/(2*(i+nf)/x-t);
a = t;
b = 1.0f;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = nf*logf(fabsf(w));
if (tmp < 88.721679688f) {
for (i=nm1; i>0; i--) {
temp = b;
b = 2.0f*i*b/x - a;
a = temp;
}
} else {
for (i=nm1; i>0; i--){
temp = b;
b = 2.0f*i*b/x - a;
a = temp;
/* scale b to avoid spurious overflow */
if (b > 0x1p60f) {
a /= b;
t /= b;
b = 1.0f;
}
}
}
z = j0f(x);
w = j1f(x);
if (fabsf(z) >= fabsf(w))
b = t*z/b;
else
b = t*w/a;
}
}
return sign ? -b : b;
}
float ynf(int n, float x)
{
uint32_t ix, ib;
int nm1, sign, i;
float a, b, temp;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix > 0x7f800000) /* nan */
return x;
if (sign && ix != 0) /* x < 0 */
return 0/0.0f;
if (ix == 0x7f800000)
return 0.0f;
if (n == 0)
return y0f(x);
if (n < 0) {
nm1 = -(n+1);
sign = n&1;
} else {
nm1 = n-1;
sign = 0;
}
if (nm1 == 0)
return sign ? -y1f(x) : y1f(x);
a = y0f(x);
b = y1f(x);
/* quit if b is -inf */
GET_FLOAT_WORD(ib,b);
for (i = 0; i < nm1 && ib != 0xff800000; ) {
i++;
temp = b;
b = (2.0f*i/x)*b - a;
GET_FLOAT_WORD(ib, b);
a = temp;
}
return sign ? -b : b;
}

View file

@ -3,8 +3,6 @@
#if !(__ASSEMBLER__ + __LINKER__ + 0)
COSMOPOLITAN_C_START_
extern int __signgam;
double __cos(double, double) hidden;
double __sin(double, double, int) hidden;
double __tan(double, double, int) hidden;

View file

@ -23,5 +23,5 @@
* Returns natural logarithm of absolute value of gamma function.
*/
double lgamma(double x) {
return lgamma_r(x, &__signgam);
return lgamma_r(x, &signgam);
}

View file

@ -20,5 +20,5 @@
#include "libc/tinymath/kernel.internal.h"
float lgammaf(float x) {
return lgammaf_r(x, &__signgam);
return lgammaf_r(x, &signgam);
}

View file

@ -37,7 +37,7 @@ asm(".ident\t\"\\n\\n\
Double-precision math functions (MIT License)\\n\
Copyright 2018 ARM Limited\"");
asm(".include \"libc/disclaimer.inc\"");
/* clang-format off */
// clang-format off
/*
* Copyright (c) 2017-2018, Arm Limited.

View file

@ -1,67 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/nexgen32e/x86feature.h"
#include "libc/intrin/smmintrin.internal.h"
#include "libc/macros.internal.h"
// Rounds to nearest integer.
//
// @param is double passed in %xmm0
// @return double in %xmm0
// @note rounding behavior can be changed in mxcsr
rint:
#if !X86_NEED(SSE4_2)
testb X86_HAVE(SSE4_2)+kCpuids(%rip)
jz rint$k8
.text.antiquity
rint$k8:
0: movq %xmm0,%rax
movq %xmm0,%rdx
shr $52,%rdx
and $2047,%edx
cmp $1074,%edx
jg 2f
movsd mmm(%rip),%xmm1
shr $63,%rax
jne 3f
addsd %xmm1,%xmm0
subsd %xmm1,%xmm0
1: pxor %xmm2,%xmm2
ucomisd %xmm2,%xmm0
jp 2f
jne 2f
movsd sgn(%rip),%xmm0
test %rax,%rax
je 4f
2: ret
3: subsd %xmm1,%xmm0
addsd %xmm1,%xmm0
jmp 1b
4: pxor %xmm0,%xmm0
ret
.endfn rint$k8,globl,hidden
.previous
.rodata.cst8
sgn: .quad 0x8000000000000000
mmm: .quad 0x4330000000000000
.previous
#endif
roundsd $_MM_FROUND_RINT,%xmm0,%xmm0
ret
.endfn rint,globl

59
libc/tinymath/rint.c Normal file
View file

@ -0,0 +1,59 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
#define EPS DBL_EPSILON
#elif FLT_EVAL_METHOD==2
#define EPS LDBL_EPSILON
#endif
static const double_t toint = 1/EPS;
double rint(double x)
{
union {double f; uint64_t i;} u = {x};
int e = u.i>>52 & 0x7ff;
int s = u.i>>63;
double_t y;
if (e >= 0x3ff+52)
return x;
if (s)
y = x - toint + toint;
else
y = x + toint - toint;
if (y == 0)
return s ? -0.0 : 0;
return y;
}

View file

@ -1,47 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
rintf: .leafprologue
.profilable
movaps %xmm0,%xmm1
movss .LC8(%rip),%xmm2
andps %xmm2,%xmm1
movss .LC7(%rip),%xmm3
comiss %xmm1,%xmm3
jbe 1f
addss %xmm3,%xmm1
andnps %xmm0,%xmm2
movaps %xmm2,%xmm0
subss %xmm3,%xmm1
orps %xmm1,%xmm0
1: .leafepilogue
.endfn rintf,globl
.rodata.cst4
.LC7: .long 1258291200
.rodata.cst16
.LC8: .long 2147483647
.long 0
.long 0
.long 0
// TODO(jart):
// vroundss $4,%xmm0,%xmm0,%xmm0

61
libc/tinymath/rintf.c Normal file
View file

@ -0,0 +1,61 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD==0
#define EPS FLT_EPSILON
#elif FLT_EVAL_METHOD==1
#define EPS DBL_EPSILON
#elif FLT_EVAL_METHOD==2
#define EPS LDBL_EPSILON
#endif
static const float_t toint = 1/EPS;
float rintf(float x)
{
union {float f; uint32_t i;} u = {x};
int e = u.i>>23 & 0xff;
int s = u.i>>31;
float_t y;
if (e >= 0x7f+23)
return x;
if (s)
y = x - toint + toint;
else
y = x + toint - toint;
if (y == 0)
return s ? -0.0f : 0.0f;
return y;
}

View file

@ -1,68 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
#include "libc/intrin/smmintrin.internal.h"
#include "libc/nexgen32e/x86feature.h"
// Rounds to nearest integer, away from zero.
//
// @param 𝑥 is double scalar in low half of %xmm0
// @return double scalar in low half of %xmm0
// @define round(𝑥) = copysign(trunc(fabs(𝑥)+.5),𝑥)
// round(𝑥) = trunc(𝑥+copysign(.5,𝑥))
round:
#if !X86_NEED(SSE4_2)
testb X86_HAVE(SSE4_2)+kCpuids(%rip)
jz round$k8
.text.antiquity
round$k8:
.leafprologue
.profilable
movapd %xmm0,%xmm1
movsd D(%rip),%xmm2
movsd C(%rip),%xmm3
andpd %xmm2,%xmm1
ucomisd %xmm1,%xmm3
jbe 2f
addsd A(%rip),%xmm1
andnpd %xmm0,%xmm2
movapd %xmm2,%xmm0
cvttsd2siq %xmm1,%rax
pxor %xmm1,%xmm1
cvtsi2sdq %rax,%xmm1
orpd %xmm1,%xmm0
2: .leafepilogue
.endfn round$k8,globl,hidden
.previous
.rodata.cst16
C: .quad 0x4330000000000000,0
D: .quad 0x7fffffffffffffff,0
.previous
#endif
movapd %xmm0,%xmm1
andpd B(%rip),%xmm0
orpd A(%rip),%xmm0
addsd %xmm1,%xmm0
roundsd $_MM_FROUND_TO_ZERO,%xmm0,%xmm0
ret
.endfn round,globl
.rodata.cst16
A: .quad 0x3fdfffffffffffff,0
B: .quad 0x8000000000000000,0

69
libc/tinymath/round.c Normal file
View file

@ -0,0 +1,69 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/internal.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
#define EPS DBL_EPSILON
#elif FLT_EVAL_METHOD==2
#define EPS LDBL_EPSILON
#endif
static const double_t toint = 1/EPS;
double round(double x)
{
union {double f; uint64_t i;} u = {x};
int e = u.i >> 52 & 0x7ff;
double_t y;
if (e >= 0x3ff+52)
return x;
if (u.i >> 63)
x = -x;
if (e < 0x3ff-1) {
/* raise inexact if x!=0 */
FORCE_EVAL(x + toint);
return 0*u.f;
}
y = x + toint - toint - x;
if (y > 0.5)
y = y + x - 1;
else if (y <= -0.5)
y = y + x + 1;
else
y = y + x;
if (u.i >> 63)
y = -y;
return y;
}

View file

@ -1,66 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
#include "libc/intrin/smmintrin.internal.h"
#include "libc/nexgen32e/x86feature.h"
// Rounds to nearest integer, away from zero.
//
// @param 𝑥 is float scalar in low quarter of %xmm0
// @return float scalar in low quarter of %xmm0
roundf:
#if !X86_NEED(SSE4_2)
testb X86_HAVE(SSE4_2)+kCpuids(%rip)
jz roundf$k8
.text.antiquity
roundf$k8:
.leafprologue
.profilable
movaps %xmm0,%xmm1
movss D(%rip),%xmm2
movss C(%rip),%xmm3
andps %xmm2,%xmm1
ucomiss %xmm1,%xmm3
jbe 2f
addss A(%rip),%xmm1
cvttss2sil %xmm1,%eax
pxor %xmm1,%xmm1
cvtsi2ssl %eax,%xmm1
andnps %xmm0,%xmm2
movaps %xmm2,%xmm0
orps %xmm1,%xmm0
2: .leafepilogue
.endfn roundf$k8,globl,hidden
.previous
.rodata.cst16
C: .long 0x4b000000,0,0,0
D: .long 0x7fffffff,0,0,0
.previous
#endif
movaps %xmm0,%xmm1
andps B(%rip),%xmm0
orps A(%rip),%xmm0
addss %xmm1,%xmm0
roundss $_MM_FROUND_TO_ZERO,%xmm0,%xmm0
ret
.endfn roundf,globl
.rodata.cst16
A: .long 0x3effffff,0,0,0
B: .long 0x80000000,0,0,0

70
libc/tinymath/roundf.c Normal file
View file

@ -0,0 +1,70 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/internal.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
#if FLT_EVAL_METHOD==0
#define EPS FLT_EPSILON
#elif FLT_EVAL_METHOD==1
#define EPS DBL_EPSILON
#elif FLT_EVAL_METHOD==2
#define EPS LDBL_EPSILON
#endif
static const float_t toint = 1/EPS;
float roundf(float x)
{
union {float f; uint32_t i;} u = {x};
int e = u.i >> 23 & 0xff;
float_t y;
if (e >= 0x7f+23)
return x;
if (u.i >> 31)
x = -x;
if (e < 0x7f-1) {
FORCE_EVAL(x + toint);
return 0*u.f;
}
y = x + toint - toint - x;
if (y > 0.5f)
y = y + x - 1;
else if (y <= -0.5f)
y = y + x + 1;
else
y = y + x;
if (u.i >> 31)
y = -y;
return y;
}

View file

@ -1,29 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
// Returns 𝑥 × 2ʸ.
//
// @param 𝑥 is double passed in %xmm0
// @param 𝑦 is double passed in %xmm1, which is truncated
// @return result in %xmm0
// @see ldexp()
scalb: cvttsd2si %xmm1,%edi
jmp ldexp
.endfn scalb,globl

70
libc/tinymath/scalb.c Normal file
View file

@ -0,0 +1,70 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2020 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
asm(".ident\t\"\\n\\n\
fdlibm (fdlibm license)\\n\
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/e_scalb.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* scalb(x, fn) is provide for
* passing various standard test suite. One
* should use scalbn() instead.
*/
double scalb(double x, double fn)
{
if (isnan(x) || isnan(fn))
return x*fn;
if (!isfinite(fn)) {
if (fn > 0.0)
return x*fn;
else
return x/(-fn);
}
if (rint(fn) != fn) return (fn-fn)/(fn-fn);
if ( fn > 65000.0) return scalbn(x, 65000);
if (-fn > 65000.0) return scalbn(x,-65000);
return scalbn(x,(int)fn);
}

View file

@ -17,4 +17,4 @@
PERFORMANCE OF THIS SOFTWARE.
*/
int __signgam;
int signgam;

View file

@ -1,28 +0,0 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi
Copyright 2020 Justine Alexandra Roberts Tunney
Permission to use, copy, modify, and/or distribute this software for
any purpose with or without fee is hereby granted, provided that the
above copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/
#include "libc/macros.internal.h"
// Returns tangent of 𝑥.
//
// @param 𝑥 is float scalar in low quarter of %xmm0
// @return float scalar in low quarter of %xmm0
// @domain -(3π/8) < 𝑥 < 3π/8 for best accuracy
tanf: ezlea tanl,ax
jmp _f2ld2
.endfn tanf,globl

99
libc/tinymath/tanf.c Normal file
View file

@ -0,0 +1,99 @@
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│
vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi
Musl Libc
Copyright © 2005-2020 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "libc/math.h"
#include "libc/tinymath/complex.internal.h"
#include "libc/tinymath/kernel.internal.h"
asm(".ident\t\"\\n\\n\
Musl libc (MIT License)\\n\
Copyright 2005-2014 Rich Felker, et. al.\"");
asm(".include \"libc/disclaimer.inc\"");
// clang-format off
/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Small multiples of pi/2 rounded to double precision. */
static const double
t1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
t2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
t3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
t4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
float tanf(float x)
{
double y;
uint32_t ix;
unsigned n, sign;
GET_FLOAT_WORD(ix, x);
sign = ix >> 31;
ix &= 0x7fffffff;
if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact if x!=0 and underflow if subnormal */
FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f);
return x;
}
return __tandf(x, 0);
}
if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */
if (ix <= 0x4016cbe3) /* |x| ~<= 3pi/4 */
return __tandf((sign ? x+t1pio2 : x-t1pio2), 1);
else
return __tandf((sign ? x+t2pio2 : x-t2pio2), 0);
}
if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */
if (ix <= 0x40afeddf) /* |x| ~<= 7*pi/4 */
return __tandf((sign ? x+t3pio2 : x-t3pio2), 1);
else
return __tandf((sign ? x+t4pio2 : x-t4pio2), 0);
}
/* tan(Inf or NaN) is NaN */
if (ix >= 0x7f800000)
return x - x;
/* argument reduction */
n = __rem_pio2f(x, &y);
return __tandf(y, n&1);
}