Attachment #54631 for bug #82654




	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
	k_cos.c k_cosf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c k_sinf.c \
	k_tan.c k_tanf.c \
	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c s_cbrtl.c \
	s_ceil.c s_ceilf.c s_ceill.c \
	s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \
	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \

	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_logl.c s_log10l.c \
        s_lrint.c s_lrintf.c \
	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
	s_nearbyint.c s_nextafter.c s_nextafterf.c \
	s_nexttowardf.c s_remquo.c s_remquof.c \
	s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
        s_sqrtl.c s_tan.c \
	s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c s_truncl.c \
	w_cabs.c w_cabsf.c w_drem.c w_dremf.c


Lines 45-52 Link Here

(-)msun/man/exp.3 (-6 / +14 lines)
45	.Nm expm1f ,	45	.Nm expm1f ,
46	.Nm log ,	46	.Nm log ,
47	.Nm logf ,	47	.Nm logf ,
		48	.Nm logl ,
48	.Nm log10 ,	49	.Nm log10 ,
49	.Nm log10f ,	50	.Nm log10f ,
		51	.Nm log10l ,
50	.Nm log1p ,	52	.Nm log1p ,
51	.Nm log1pf ,	53	.Nm log1pf ,
52	.Nm pow ,	54	.Nm pow ,
Lines 72-81 Link Here
72	.Fn log "double x"	74	.Fn log "double x"
73	.Ft float	75	.Ft float
74	.Fn logf "float x"	76	.Fn logf "float x"
		77	.Ft "long double"
		78	.Fn logl "long double x"
75	.Ft double	79	.Ft double
76	.Fn log10 "double x"	80	.Fn log10 "double x"
77	.Ft float	81	.Ft float
78	.Fn log10f "float x"	82	.Fn log10f "float x"
		83	.Ft "long double"
		84	.Fn log10l "long double x"
79	.Ft double	85	.Ft double
80	.Fn log1p "double x"	86	.Fn log1p "double x"
81	.Ft float	87	.Ft float
Lines 109-124 Link Here
109	.Fa x .	115	.Fa x .
110	.Pp	116	.Pp
111	The	117	The
112	.Fn log	118	.Fn log ,
113	and the	119	.Fn logf ,
114	.Fn logf	120	and
		121	.Fn logl
115	functions compute the value of the natural logarithm of argument	122	functions compute the value of the natural logarithm of argument
116	.Fa x .	123	.Fa x .
117	.Pp	124	.Pp
118	The	125	The
119	.Fn log10	126	.Fn log10 ,
120	and the	127	.Fn log10f ,
121	.Fn log10f	128	and
		129	.Fn log10l
122	functions compute the value of the logarithm of argument	130	functions compute the value of the logarithm of argument
123	.Fa x	131	.Fa x
124	to base 10.	132	to base 10.

Lines 49-83 Link Here

(-)msun/man/sqrt.3 (-14 / +22 lines)
49	.Fn cbrt "double x"	49	.Fn cbrt "double x"
50	.Ft float	50	.Ft float
51	.Fn cbrtf "float x"	51	.Fn cbrtf "float x"
		52	.Ft "long double"
		53	.Fn cbrtl "long double x"
52	.Ft double	54	.Ft double
53	.Fn sqrt "double x"	55	.Fn sqrt "double x"
54	.Ft float	56	.Ft float
55	.Fn sqrtf "float x"	57	.Fn sqrtf "float x"
		58	.Ft "long double"
		59	.Fn sqrtl "long double x"
56	.Sh DESCRIPTION	60	.Sh DESCRIPTION
57	The	61	The
58	.Fn cbrt	62	.Fn cbrt ,
59	and the	63	.Fn cbrtf ,
60	.Fn cbrtf	64	and
		65	.Fn cbrtl
61	functions compute	66	functions compute
62	the cube root of	67	the cube root of
63	.Ar x .	68	.Ar x .
64	.Pp	69	.Pp
65	The	70	The
66	.Fn sqrt	71	.Fn sqrt ,
67	and the	72	.Fn sqrtf ,
68	.Fn sqrtf	73	and
69	functions compute the	74	.Fn sqrtl
70	non-negative square root of x.	75	functions compute the non-negative square root of
		76	.Ar x .
71	.Sh RETURN VALUES	77	.Sh RETURN VALUES
72	The	78	The
73	.Fn cbrt	79	.Fn cbrt ,
74	and the	80	.Fn cbrtf ,
75	.Fn cbrtf	81	and
		82	.Fn cbrtl
76	functions return the requested cube root.	83	functions return the requested cube root.
77	The	84	The
78	.Fn sqrt	85	.Fn sqrt ,
79	and the	86	.Fn sqrtf ,
80	.Fn sqrtf	87	and
		88	.Fn sqrtl
81	functions return the requested square root	89	functions return the requested square root
82	unless an error occurs.	90	unless an error occurs.
83	An attempt to take the	91	An attempt to take the




long double	atan2l(long double, long double);
long double	atanhl(long double);
long double	atanl(long double);
long double	cbrtl(long double);
#endif
long double	cbrtl(long double);
long double	ceill(long double);
long double	copysignl(long double, long double) __pure2;
#if 0

long long	llrintl(long double);
#endif
long long	llroundl(long double);
#if 0
long double	log10l(long double);
#if 0
long double	log1pl(long double);
long double	log2l(long double);
long double	logbl(long double);
#endif
long double	logl(long double);
#if 0
long		lrintl(long double);
#endif
long		lroundl(long double);

#if 0
long double	sinhl(long double);
long double	sinl(long double);
#endif
long double	sqrtl(long double);
#if 0
long double	tanhl(long double);
long double	tanl(long double);
long double	tgammal(long double);




/*-
 * Copyright (c) 2005, Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <math.h>
#include <float.h>

/*
 *  Compute the cube root of a long double argument by decomposing the
 *  the argument into its base 2 fraction and exponent.
 *
 *  x**(1/3) = (f 2**n)**(1/3)
 *           = f**(1/3) * 2**(n/3)                 (mod(n,3) = 0)
 *           = f**(1/3) * 2**(1/3) * 2**[(n-1)/3]  (mod(n,3) = 1)
 *           = f**(1/3) * 2**(2/3) * 2**[(n-2)/3]  (mod(n,3) = 2)
 *
 *  where f = [1/2,1).
 *
 *  Use Newton's rule to compute f**(1/3)  
 *     y_(k+1) = (x / y_k**2 + 2*y_k) / 3 
 *  where k is the iteration count. 
 */

long double cbrtl(long double x) {

	int i, n, s = 1;
	long double f, yk, oyk;

	/* cbrtl(x) = NaN or +-Inf. */
	if (isinf(x) || isnan(x))
		return (x+x);

	/* Save the sign. */
	if (x < 0.L) {
		s = -1;
		x = -x;
	}

	f = frexpl(x, &n);

	yk = oyk = f;
	while(1) {
		yk = (f / yk / yk + 2.0L * yk) / 3.0L;
		if (fabsl(oyk - yk) < LDBL_EPSILON) break;
		oyk = yk;
	}

	switch (n%3) {
		case 0:
			yk = ldexpl(yk, n / 3);
			break;
		case 1:
			yk = ldexpl(yk, (n-1) / 3) * 1.259921049894873165L;
			break;
		case 2:
			yk = ldexpl(yk, (n-2) / 3) * 1.587401051968199475L;
			break;
	}

	return (s * yk);

}




/*-
 * Copyright (c) 2005, Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 *  Compute the base 10 logrithm of the argument x by .
 *
 *  log10(x) = log(x) / log(10).
 *
 *  where log(x) is computed via logl(x).
 */

#include <math.h>

#define LOG10L2 3.010299956639811952e-1L
#define LOGL10	2.302585092994045684L

static long double zero = 0.0L;

long double log10l(long double x) {

	int n;
	long double f, val;

	/* log10(x) = sNaN if x < 0. */
	if (x < 0.0L)
		return (x - x) / (x - x);
 
	/* log10(+Inf) = +Inf */
	if (isinf(x))
		return x*x+x;

	/* log10(+-0) = -Inf with signal */
	if (x == 0.0L)
		return (- 1.0L / zero);

	/* log10(NaN) = NaN */
	if (isnan(x))
		return x;

	val = logl(x) / LOGL10;

	return val;

}




/*-
 * Copyright (c) 2005, Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 *  Compute the natural logrithm of x by decomposing x into its base 2
 *  representation.
 *
 *  log(x) = log(f * 2**n)
 *         = log(f) + n * log(2)
 *
 *  where f = [1/2,1).
 *
 * Use a Taylor's series about f0 = 0.75 to compute log(f).
 *
 *  log(f) = log(f0) + sum{(1/n) * (-1)**(n-1) * ((f - f0)/f0)**n}.
 */

#include <math.h>
#include <float.h>

#define LOGL2	 6.931471805599453094e-1L
#define LOGLF0	-2.876820724517809274e-1L

static long double zero = 0.0L;

long double logl(long double x) {

	int i, n;
	long double f, t, val, oval;

	/* log(x) = sNaN if x < 0. */
	if (x < 0.0L)
		return (x - x) / (x - x);
 
	/* log(+Inf) = +Inf */
	if (isinf(x))
		return x*x+x;

	/* log(+-0) = -Inf with signal */
	if (x == 0.0L)
		return (- 1.0L / zero);

	/* log(NaN) =  NaN with signal */
	if (isnan(x))
		return (x+x);

	/*
	 *  Special case for values near log(1).  Use the series expansion
	 *  for log(1+x) = x - (1/2) * x**2 + (1/3) * x**3 + ...
	 */
	if (0.95L < x && x < 1.05L) {

		f = t = x - 1.0L;
		oval = val = 0.L;

		for (i = 1; ; i++) {
			val += t / (long double) i;
			t *= -f;
			if (fabsl(oval - val) < LDBL_EPSILON) break;
			oval = val;
		} 

		return (val);
	}

	f = frexpl(x, &n);
	f = t = (f - 0.75L) / 0.75L;

	oval = val = LOGLF0;
	for (i = 1; ; i++) {
		val += t / (long double) i;
		t *= -f;
		if (fabsl(oval-val) < LDBL_EPSILON) break;
		oval = val;
	} 

	return (val + n * LOGL2);

}




/*-
 * Copyright (c) 2005, Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 *  Compute the sqrt root of a long double argument by decomposing the
 *  the argument into its base 2 fraction and exponent.
 *
 *  x**(1/2) = (f 2**n)**(1/2)
 *           = f**(1/2) * 2**(n/2)                 (n even)
 *           = f**(1/2) * 2**(1/2) * 2**[(n-1)/2]  (n odd)
 *  where f = [1/2,1).
 *
 *  Use Newton's rule to compute f
 *     y_(k+1) = (1/2) * (x / y_k + y_k)
 *  where k is the iteration count.
 */

#include <math.h>
#include <float.h>

long double sqrtl(long double x) {

	int n;
	long double f, yk, oyk;

	/* sqrt(NaN) = NaN, sqrt(+Inf) = +Inf, or  sqrt(-Inf) = sNaN */
	if (isnan(x) || isinf(x))
		return x*x+x;
		
	/* sqrt(+-0) = 0 */
	if (x == 0.0L)
		return fabsl(x);

	/* sqrt(x), x < 0 */
	if (x < 0.0L)
		return (x - x) / (x - x);

	f = frexpl(x, &n);

	yk = oyk = f;
	while(1) {
		yk = 0.5L * (f / yk + yk);
		if (fabsl(oyk - yk) < LDBL_EPSILON) break;
		oyk = yk;
	}

	if (n%2 == 0) 	
		yk = ldexpl(yk, n / 2);
	else
		yk = 1.41421356237309505L * ldexpl(yk,(n - 1) / 2);

	return yk;

}

Return to bug 82654

Lines 38-44 Link Here

(-)msun/Makefile (-3 / +5 lines)
38	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \	38	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
39	k_cos.c k_cosf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c k_sinf.c \	39	k_cos.c k_cosf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c k_sinf.c \
40	k_tan.c k_tanf.c \	40	k_tan.c k_tanf.c \
41	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c \	41	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c s_cbrtl.c \
42	s_ceil.c s_ceilf.c s_ceill.c \	42	s_ceil.c s_ceilf.c s_ceill.c \
43	s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \	43	s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \
44	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \	44	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
Lines 48-60 Link Here
48	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \	48	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
49	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \	49	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
50	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \	50	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
51	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \	51	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_logl.c s_log10l.c \
		52	s_lrint.c s_lrintf.c \
52	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \	53	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
53	s_nearbyint.c s_nextafter.c s_nextafterf.c \	54	s_nearbyint.c s_nextafter.c s_nextafterf.c \
54	s_nexttowardf.c s_remquo.c s_remquof.c \	55	s_nexttowardf.c s_remquo.c s_remquof.c \
55	s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \	56	s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
56	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \	57	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
57	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c s_tan.c \	58	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
		59	s_sqrtl.c s_tan.c \
58	s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c s_truncl.c \	60	s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c s_truncl.c \
59	w_cabs.c w_cabsf.c w_drem.c w_dremf.c	61	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
60		62

Lines 397-404 Link Here

(-)msun/src/math.h (-2 / +6 lines)
397	long double atan2l(long double, long double);	397	long double atan2l(long double, long double);
398	long double atanhl(long double);	398	long double atanhl(long double);
399	long double atanl(long double);	399	long double atanl(long double);
400	long double cbrtl(long double);
401	#endif	400	#endif
		401	long double cbrtl(long double);
402	long double ceill(long double);	402	long double ceill(long double);
403	long double copysignl(long double, long double) __pure2;	403	long double copysignl(long double, long double) __pure2;
404	#if 0	404	#if 0
Lines 430-441 Link Here
430	long long llrintl(long double);	430	long long llrintl(long double);
431	#endif	431	#endif
432	long long llroundl(long double);	432	long long llroundl(long double);
433	#if 0
434	long double log10l(long double);	433	long double log10l(long double);
		434	#if 0
435	long double log1pl(long double);	435	long double log1pl(long double);
436	long double log2l(long double);	436	long double log2l(long double);
437	long double logbl(long double);	437	long double logbl(long double);
		438	#endif
438	long double logl(long double);	439	long double logl(long double);
		440	#if 0
439	long lrintl(long double);	441	long lrintl(long double);
440	#endif	442	#endif
441	long lroundl(long double);	443	long lroundl(long double);
Lines 460-466 Link Here
460	#if 0	462	#if 0
461	long double sinhl(long double);	463	long double sinhl(long double);
462	long double sinl(long double);	464	long double sinl(long double);
		465	#endif
463	long double sqrtl(long double);	466	long double sqrtl(long double);
		467	#if 0
464	long double tanhl(long double);	468	long double tanhl(long double);
465	long double tanl(long double);	469	long double tanl(long double);
466	long double tgammal(long double);	470	long double tgammal(long double);

Line 0 Link Here

(-)msun/src/s_cbrtl.c (+84 lines)
		1	/*-
		2	* Copyright (c) 2005, Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	#include <math.h>
		28	#include <float.h>
		29
		30	/*
		31	* Compute the cube root of a long double argument by decomposing the
		32	* the argument into its base 2 fraction and exponent.
		33	*
		34	* x(1/3) = (f 2n)**(1/3)
		35	* = f*(1/3) 2**(n/3) (mod(n,3) = 0)
		36	* = f*(1/3) 2*(1/3) 2**[(n-1)/3] (mod(n,3) = 1)
		37	* = f*(1/3) 2*(2/3) 2**[(n-2)/3] (mod(n,3) = 2)
		38	*
		39	* where f = [1/2,1).
		40	*
		41	* Use Newton's rule to compute f**(1/3)
		42	* y_(k+1) = (x / y_k*2 + 2y_k) / 3
		43	* where k is the iteration count.
		44	*/
		45
		46	long double cbrtl(long double x) {
		47
		48	int i, n, s = 1;
		49	long double f, yk, oyk;
		50
		51	/* cbrtl(x) = NaN or +-Inf. */
		52	if (isinf(x) \|\| isnan(x))
		53	return (x+x);
		54
		55	/* Save the sign. */
		56	if (x < 0.L) {
		57	s = -1;
		58	x = -x;
		59	}
		60
		61	f = frexpl(x, &n);
		62
		63	yk = oyk = f;
		64	while(1) {
65	yk = (f / yk / yk + 2.0L * yk) / 3.0L;
66	if (fabsl(oyk - yk) < LDBL_EPSILON) break;
67	oyk = yk;
68	}
69
70	switch (n%3) {
71	case 0:
72	yk = ldexpl(yk, n / 3);
73	break;
74	case 1:
75	yk = ldexpl(yk, (n-1) / 3) * 1.259921049894873165L;
76	break;
77	case 2:
78	yk = ldexpl(yk, (n-2) / 3) * 1.587401051968199475L;
79	break;
80	}
81
82	return (s * yk);
83
84	}

Line 0 Link Here

(-)msun/src/s_log10l.c (+67 lines)
		1	/*-
		2	* Copyright (c) 2005, Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	/*
		28	* Compute the base 10 logrithm of the argument x by .
		29	*
		30	* log10(x) = log(x) / log(10).
		31	*
		32	* where log(x) is computed via logl(x).
		33	*/
		34
		35	#include <math.h>
		36
		37	#define LOG10L2 3.010299956639811952e-1L
		38	#define LOGL10 2.302585092994045684L
		39
		40	static long double zero = 0.0L;
		41
		42	long double log10l(long double x) {
		43
		44	int n;
		45	long double f, val;
		46
		47	/* log10(x) = sNaN if x < 0. */
		48	if (x < 0.0L)
		49	return (x - x) / (x - x);
		50
		51	/* log10(+Inf) = +Inf */
		52	if (isinf(x))
		53	return x*x+x;
		54
		55	/* log10(+-0) = -Inf with signal */
		56	if (x == 0.0L)
		57	return (- 1.0L / zero);
		58
		59	/* log10(NaN) = NaN */
		60	if (isnan(x))
		61	return x;
		62
		63	val = logl(x) / LOGL10;
		64
65	return val;
66
67	}

Line 0 Link Here

(-)msun/src/s_logl.c (+102 lines)
		1	/*-
		2	* Copyright (c) 2005, Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	/*
		28	* Compute the natural logrithm of x by decomposing x into its base 2
		29	* representation.
		30	*
		31	* log(x) = log(f * 2**n)
		32	* = log(f) + n * log(2)
		33	*
		34	* where f = [1/2,1).
		35	*
		36	* Use a Taylor's series about f0 = 0.75 to compute log(f).
		37	*
		38	* log(f) = log(f0) + sum{(1/n) * (-1)*(n-1) ((f - f0)/f0)**n}.
		39	*/
		40
		41	#include <math.h>
		42	#include <float.h>
		43
		44	#define LOGL2 6.931471805599453094e-1L
		45	#define LOGLF0 -2.876820724517809274e-1L
		46
		47	static long double zero = 0.0L;
		48
		49	long double logl(long double x) {
		50
		51	int i, n;
		52	long double f, t, val, oval;
		53
		54	/* log(x) = sNaN if x < 0. */
		55	if (x < 0.0L)
		56	return (x - x) / (x - x);
		57
		58	/* log(+Inf) = +Inf */
		59	if (isinf(x))
		60	return x*x+x;
		61
		62	/* log(+-0) = -Inf with signal */
		63	if (x == 0.0L)
		64	return (- 1.0L / zero);
65
66	/* log(NaN) = NaN with signal */
67	if (isnan(x))
68	return (x+x);
69
70	/*
71	* Special case for values near log(1). Use the series expansion
72	* for log(1+x) = x - (1/2) * x*2 + (1/3) x**3 + ...
73	*/
74	if (0.95L < x && x < 1.05L) {
75
76	f = t = x - 1.0L;
77	oval = val = 0.L;
78
79	for (i = 1; ; i++) {
80	val += t / (long double) i;
81	t *= -f;
82	if (fabsl(oval - val) < LDBL_EPSILON) break;
83	oval = val;
84	}
85
86	return (val);
87	}
88
89	f = frexpl(x, &n);
90	f = t = (f - 0.75L) / 0.75L;
91
92	oval = val = LOGLF0;
93	for (i = 1; ; i++) {
94	val += t / (long double) i;
95	t *= -f;
96	if (fabsl(oval-val) < LDBL_EPSILON) break;
97	oval = val;
98	}
99
100	return (val + n * LOGL2);
101
102	}

Line 0 Link Here

(-)msun/src/s_sqrtl.c (+77 lines)
		1	/*-
		2	* Copyright (c) 2005, Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	/*
		28	* Compute the sqrt root of a long double argument by decomposing the
		29	* the argument into its base 2 fraction and exponent.
		30	*
		31	* x(1/2) = (f 2n)**(1/2)
		32	* = f*(1/2) 2**(n/2) (n even)
		33	* = f*(1/2) 2*(1/2) 2**[(n-1)/2] (n odd)
		34	* where f = [1/2,1).
		35	*
		36	* Use Newton's rule to compute f
		37	* y_(k+1) = (1/2) * (x / y_k + y_k)
		38	* where k is the iteration count.
		39	*/
		40
		41	#include <math.h>
		42	#include <float.h>
		43
		44	long double sqrtl(long double x) {
		45
		46	int n;
		47	long double f, yk, oyk;
		48
		49	/* sqrt(NaN) = NaN, sqrt(+Inf) = +Inf, or sqrt(-Inf) = sNaN */
		50	if (isnan(x) \|\| isinf(x))
		51	return x*x+x;
		52
		53	/* sqrt(+-0) = 0 */
		54	if (x == 0.0L)
		55	return fabsl(x);
		56
		57	/* sqrt(x), x < 0 */
		58	if (x < 0.0L)
		59	return (x - x) / (x - x);
		60
		61	f = frexpl(x, &n);
		62
		63	yk = oyk = f;
		64	while(1) {
65	yk = 0.5L * (f / yk + yk);
66	if (fabsl(oyk - yk) < LDBL_EPSILON) break;
67	oyk = yk;
68	}
69
70	if (n%2 == 0)
71	yk = ldexpl(yk, n / 2);
72	else
73	yk = 1.41421356237309505L * ldexpl(yk,(n - 1) / 2);
74
75	return yk;
76
77	}