Update to FreeBSD libm r336665.
This reverts commit 253a8306316cedfd6fd3e3a169fbffe4cac04035 and moves
us forward to a revision that contains fixes for the problem with the
previous attempt.
This also makes sincos(3)/sincosf(3)/sincosl(3) available to `_BSD_SOURCE`
as well as `_GNU_SOURCE`.
The new FreeBSD libm code requires the FreeBSD `__CONCAT` macro, and all
our existing callers are FreeBSD too, so update that.
There's also an assumption that <complex.h> drags in <math.h> which isn't
true for us, so work around that with `-include` in the makefile. This
then causes clang to recognize a bug -- returning from a void function --
in our fake (LP32) sincosl(3), so fix that too.
Bug: http://b/111710419
Change-Id: I84703ad844f8afde6ec6b11604ab3c096ccb62c3
Test: ran tests
diff --git a/libm/upstream-freebsd/lib/msun/src/e_jn.c b/libm/upstream-freebsd/lib/msun/src/e_jn.c
index a1130c5..58ec905 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_jn.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_jn.c
@@ -1,4 +1,3 @@
-
/* @(#)e_jn.c 1.4 95/01/18 */
/*
* ====================================================
@@ -6,19 +5,19 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_jn.c 279856 2015-03-10 17:10:54Z kargl $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_jn.c 336089 2018-07-08 16:26:13Z markj $");
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
- *
+ *
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
@@ -37,7 +36,6 @@
* yn(n,x) is similar in all respects, except
* that forward recursion is used for all
* values of n>1.
- *
*/
#include "math.h"
@@ -66,7 +64,7 @@
ix = 0x7fffffff&hx;
/* if J(n,NaN) is NaN */
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
- if(n<0){
+ if(n<0){
n = -n;
x = -x;
hx ^= 0x80000000;
@@ -77,14 +75,14 @@
x = fabs(x);
if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
b = zero;
- else if((double)n<=x) {
+ else if((double)n<=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)
+ /* (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
+ * Let s=sin(x), c=cos(x),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
@@ -100,7 +98,7 @@
case 3: temp = cos(x)-sin(x); break;
}
b = invsqrtpi*temp/sqrt(x);
- } else {
+ } else {
a = __ieee754_j0(x);
b = __ieee754_j1(x);
for(i=1;i<n;i++){
@@ -111,7 +109,7 @@
}
} else {
if(ix<0x3e100000) { /* x < 2**-29 */
- /* x is tiny, return the first Taylor expansion of J(n,x)
+ /* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if(n>33) /* underflow */
@@ -126,14 +124,14 @@
}
} else {
/* use backward recurrence */
- /* x x^2 x^2
+ /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
- * 1 1 1
+ * 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
@@ -149,9 +147,9 @@
* 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
+ * 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,v;
@@ -237,11 +235,11 @@
if(n==1) return(sign*__ieee754_y1(x));
if(ix==0x7ff00000) return zero;
if(ix>=0x52D00000) { /* x > 2**302 */
- /* (x >> n**2)
+ /* (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
+ * Let s=sin(x), c=cos(x),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------