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
      *		----------------------------------