1/*-
2 * Copyright (c) 2008-2010 David Schultz <das@FreeBSD.org>
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, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 *    notice, this list of conditions and the following disclaimer in the
12 *    documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27/*
28 * Tests for corner cases in log*().
29 */
30
31#include <sys/cdefs.h>
32__FBSDID("$FreeBSD$");
33
34#include <sys/param.h>
35#include <assert.h>
36#include <fenv.h>
37#include <float.h>
38#include <math.h>
39#include <stdio.h>
40
41#ifdef __i386__
42#include <ieeefp.h>
43#endif
44
45#include "test-utils.h"
46
47#pragma STDC FENV_ACCESS ON
48
49/*
50 * Test that a function returns the correct value and sets the
51 * exception flags correctly. The exceptmask specifies which
52 * exceptions we should check. We need to be lenient for several
53 * reasoons, but mainly because on some architectures it's impossible
54 * to raise FE_OVERFLOW without raising FE_INEXACT.
55 *
56 * These are macros instead of functions so that assert provides more
57 * meaningful error messages.
58 *
59 * XXX The volatile here is to avoid gcc's bogus constant folding and work
60 *     around the lack of support for the FENV_ACCESS pragma.
61 */
62#define	test(func, x, result, exceptmask, excepts)	do {		\
63	volatile long double _d = x;					\
64	assert(feclearexcept(FE_ALL_EXCEPT) == 0);			\
65	assert(fpequal((func)(_d), (result)));				 \
66	assert(((void)(func), fetestexcept(exceptmask) == (excepts)));	\
67} while (0)
68
69#define	test(func, x, result, exceptmask, excepts)	do {		\
70	volatile long double _d = x;					\
71	assert(feclearexcept(FE_ALL_EXCEPT) == 0);			\
72	assert(fpequal((func)(_d), (result)));				 \
73	assert(((void)(func), fetestexcept(exceptmask) == (excepts)));	\
74} while (0)
75
76#define	test_tol(func, z, result, tol)			do {		\
77	volatile long double _d = z;					\
78	debug("  testing %6s(%15La) ~= % .36Le\n", #func, _d, result);	\
79	assert(fpequal_tol((func)(_d), (result), (tol), CS_BOTH));	\
80} while (0)
81
82/* Test all the functions that compute log(x). */
83#define	testall0(x, result, exceptmask, excepts)	do {		\
84	test(log, x, result, exceptmask, excepts);			\
85	test(logf, x, result, exceptmask, excepts);			\
86	test(logl, x, result, exceptmask, excepts);			\
87	test(log2, x, result, exceptmask, excepts);			\
88	test(log2f, x, result, exceptmask, excepts);			\
89	test(log2l, x, result, exceptmask, excepts);			\
90	test(log10, x, result, exceptmask, excepts);			\
91	test(log10f, x, result, exceptmask, excepts);			\
92	test(log10l, x, result, exceptmask, excepts);			\
93} while (0)
94
95/* Test all the functions that compute log(1+x). */
96#define	testall1(x, result, exceptmask, excepts)	do {		\
97	test(log1p, x, result, exceptmask, excepts);			\
98	test(log1pf, x, result, exceptmask, excepts);			\
99	test(log1pl, x, result, exceptmask, excepts);			\
100} while (0)
101
102static void
103run_generic_tests(void)
104{
105
106	/* log(1) == 0, no exceptions raised */
107	testall0(1.0, 0.0, ALL_STD_EXCEPT, 0);
108	testall1(0.0, 0.0, ALL_STD_EXCEPT, 0);
109	testall1(-0.0, -0.0, ALL_STD_EXCEPT, 0);
110
111	/* log(NaN) == NaN, no exceptions raised */
112	testall0(NAN, NAN, ALL_STD_EXCEPT, 0);
113	testall1(NAN, NAN, ALL_STD_EXCEPT, 0);
114
115	/* log(Inf) == Inf, no exceptions raised */
116	testall0(INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
117	testall1(INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
118
119	/* log(x) == NaN for x < 0, invalid exception raised */
120	testall0(-INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
121	testall1(-INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
122	testall0(-1.0, NAN, ALL_STD_EXCEPT, FE_INVALID);
123	testall1(-1.5, NAN, ALL_STD_EXCEPT, FE_INVALID);
124
125	/* log(0) == -Inf, divide-by-zero exception */
126	testall0(0.0, -INFINITY, ALL_STD_EXCEPT & ~FE_INEXACT, FE_DIVBYZERO);
127	testall0(-0.0, -INFINITY, ALL_STD_EXCEPT & ~FE_INEXACT, FE_DIVBYZERO);
128	testall1(-1.0, -INFINITY, ALL_STD_EXCEPT & ~FE_INEXACT, FE_DIVBYZERO);
129}
130
131static void
132run_log2_tests(void)
133{
134	unsigned i;
135
136	/*
137	 * We should insist that log2() return exactly the correct
138	 * result and not raise an inexact exception for powers of 2.
139	 */
140	assert(feclearexcept(FE_ALL_EXCEPT) == 0);
141	for (i = FLT_MIN_EXP - FLT_MANT_DIG; i < FLT_MAX_EXP; i++) {
142		assert(log2f(ldexpf(1.0, i)) == i);
143		assert(fetestexcept(ALL_STD_EXCEPT) == 0);
144	}
145	for (i = DBL_MIN_EXP - DBL_MANT_DIG; i < DBL_MAX_EXP; i++) {
146		assert(log2(ldexp(1.0, i)) == i);
147		assert(fetestexcept(ALL_STD_EXCEPT) == 0);
148	}
149	for (i = LDBL_MIN_EXP - LDBL_MANT_DIG; i < LDBL_MAX_EXP; i++) {
150		assert(log2l(ldexpl(1.0, i)) == i);
151#if 0
152		/* XXX This test does not pass yet. */
153		assert(fetestexcept(ALL_STD_EXCEPT) == 0);
154#endif
155	}
156}
157
158static void
159run_roundingmode_tests(void)
160{
161
162	/*
163	 * Corner cases in other rounding modes.
164	 */
165	fesetround(FE_DOWNWARD);
166	/* These are still positive per IEEE 754R */
167#if 0
168	testall0(1.0, 0.0, ALL_STD_EXCEPT, 0);
169#else
170	/* logl, log2l, and log10l don't pass yet. */
171	test(log, 1.0, 0.0, ALL_STD_EXCEPT, 0);
172	test(logf, 1.0, 0.0, ALL_STD_EXCEPT, 0);
173	test(log2, 1.0, 0.0, ALL_STD_EXCEPT, 0);
174	test(log2f, 1.0, 0.0, ALL_STD_EXCEPT, 0);
175	test(log10, 1.0, 0.0, ALL_STD_EXCEPT, 0);
176	test(log10f, 1.0, 0.0, ALL_STD_EXCEPT, 0);
177#endif
178	testall1(0.0, 0.0, ALL_STD_EXCEPT, 0);
179	fesetround(FE_TOWARDZERO);
180	testall0(1.0, 0.0, ALL_STD_EXCEPT, 0);
181	testall1(0.0, 0.0, ALL_STD_EXCEPT, 0);
182
183	fesetround(FE_UPWARD);
184	testall0(1.0, 0.0, ALL_STD_EXCEPT, 0);
185	testall1(0.0, 0.0, ALL_STD_EXCEPT, 0);
186	/* log1p(-0.0) == -0.0 even when rounding upwards */
187	testall1(-0.0, -0.0, ALL_STD_EXCEPT, 0);
188
189	fesetround(FE_TONEAREST);
190}
191
192static void
193run_accuracy_tests(void)
194{
195	static const struct {
196		float x;
197		long double log2x;
198		long double logex;
199		long double log10x;
200        } tests[] = {
201		{  0x1p-120 + 0x1p-140,
202		  -1.19999998624139449158861798943319717e2L,
203		  -8.31776607135195754708796206665656732e1L,
204		  -3.61235990655024477716980559136055915e1L,
205		},
206		{  1.0 - 0x1p-20,
207		  -1.37586186296463416424364914705656460e-6L,
208		  -9.53674771153890007250243736279163253e-7L,
209		  -4.14175690642480911859354110516159131e-7L, },
210		{  1.0 + 0x1p-20,
211		   1.37586055084113820105668028340371476e-6L,
212		   9.53673861659188233908415514963336144e-7L,
213		   4.14175295653950611453333571759200697e-7L },
214		{  19.75,
215		   4.30378074817710292442728634194115348e0L,
216		   2.98315349134713087533848129856505779e0L,
217		   1.29556709996247903756734359702926363e0L },
218		{  19.75 * 0x1p100,
219		   1.043037807481771029244272863419411534e2L,
220		   7.229787154734166181706169344438271459e1L,
221		   3.139856666636059855894123306947856631e1L },
222	};
223        unsigned i;
224
225	for (i = 0; i < nitems(tests); i++) {
226		test_tol(log2, tests[i].x, tests[i].log2x, DBL_ULP());
227		test_tol(log2f, tests[i].x, tests[i].log2x, FLT_ULP());
228		test_tol(log2l, tests[i].x, tests[i].log2x, LDBL_ULP());
229		test_tol(log, tests[i].x, tests[i].logex, DBL_ULP());
230		test_tol(logf, tests[i].x, tests[i].logex, FLT_ULP());
231		test_tol(logl, tests[i].x, tests[i].logex, LDBL_ULP());
232		test_tol(log10, tests[i].x, tests[i].log10x, DBL_ULP());
233		test_tol(log10f, tests[i].x, tests[i].log10x, FLT_ULP());
234		test_tol(log10l, tests[i].x, tests[i].log10x, LDBL_ULP());
235		if (tests[i].x >= 0.5) {
236			test_tol(log1p, tests[i].x - 1, tests[i].logex,
237				 DBL_ULP());
238			test_tol(log1pf, tests[i].x - 1, tests[i].logex,
239				 FLT_ULP());
240			test_tol(log1pl, tests[i].x - 1, tests[i].logex,
241				 LDBL_ULP());
242		}
243	}
244}
245
246static void
247run_log1p_accuracy_tests(void)
248{
249
250	test_tol(log1pf, 0x0.333333p0F,
251		 1.82321546859847114303367992804596800640e-1L, FLT_ULP());
252	test_tol(log1p, 0x0.3333333333333p0,
253		 1.82321556793954589204283870982629267635e-1L, DBL_ULP());
254	test_tol(log1pl, 0x0.33333333333333332p0L,
255		 1.82321556793954626202683007050468762914e-1L, LDBL_ULP());
256
257	test_tol(log1pf, -0x0.333333p0F,
258		 -2.23143536413048672940940199918017467652e-1L, FLT_ULP());
259	test_tol(log1p, -0x0.3333333333333p0,
260		 -2.23143551314209700255143859052009022937e-1L, DBL_ULP());
261	test_tol(log1pl, -0x0.33333333333333332p0L,
262		 -2.23143551314209755752742563153765697950e-1L, LDBL_ULP());
263}
264
265int
266main(void)
267{
268
269	printf("1..5\n");
270
271	run_generic_tests();
272	printf("ok 1 - logarithm\n");
273
274	run_log2_tests();
275	printf("ok 2 - logarithm\n");
276
277	run_roundingmode_tests();
278	printf("ok 3 - logarithm\n");
279
280	run_accuracy_tests();
281	printf("ok 4 - logarithm\n");
282
283	run_log1p_accuracy_tests();
284	printf("ok 5 - logarithm\n");
285
286	return (0);
287}
288