hack4glarus » hack4glarus-2019-summer

k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c \

k_sinf.c k_tan.c k_tanf.c \

/*-

* SPDX-License-Identifier: BSD-2-Clause-FreeBSD

*

* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>

* 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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 <sys/cdefs.h>

#include <complex.h>

#include "math.h"

#include "math_private.h"

static const uint32_t k = 1799; /* constant for reduction */

static const double kln2 = 1246.97177782734161156; /* k * ln2 */

/*

* Compute exp(x), scaled to avoid spurious overflow. An exponent is

* returned separately in 'expt'.

*

* Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91

* Output: 2**1023 <= y < 2**1024

*/

static double

__frexp_exp(double x, int *expt)

{

double exp_x;

uint32_t hx;

/*

* We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to

* minimize |exp(kln2) - 2**k|. We also scale the exponent of

* exp_x to MAX_EXP so that the result can be multiplied by

* a tiny number without losing accuracy due to denormalization.

*/

exp_x = exp(x - kln2);

GET_HIGH_WORD(hx, exp_x);

*expt = (hx >> 20) - (0x3ff + 1023) + k;

SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));

return (exp_x);

}

/*

* __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.

* They are intended for large arguments (real part >= ln(DBL_MAX))

* where care is needed to avoid overflow.

*

* The present implementation is narrowly tailored for our hyperbolic and

* exponential functions. We assume expt is small (0 or -1), and the caller

* has filtered out very large x, for which overflow would be inevitable.

*/

double

__ldexp_exp(double x, int expt)

{

double exp_x, scale;

int ex_expt;

exp_x = __frexp_exp(x, &ex_expt);

expt += ex_expt;

INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);

return (exp_x * scale);

}

double complex

__ldexp_cexp(double complex z, int expt)

{

double x, y, exp_x, scale1, scale2;

int ex_expt, half_expt;

x = creal(z);

y = cimag(z);

exp_x = __frexp_exp(x, &ex_expt);

expt += ex_expt;

/*

* Arrange so that scale1 * scale2 == 2**expt. We use this to

* compensate for scalbn being horrendously slow.

*/

half_expt = expt / 2;

INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);

half_expt = expt - half_expt;

INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);

return (CMPLX(cos(y) * exp_x * scale1 * scale2,

sin(y) * exp_x * scale1 * scale2));

}

/*-

* SPDX-License-Identifier: BSD-2-Clause-FreeBSD

*

* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>

* 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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 <sys/cdefs.h>

#include <complex.h>

#include "math.h"

#include "math_private.h"

static const uint32_t k = 235; /* constant for reduction */

static const float kln2 = 162.88958740F; /* k * ln2 */

/*

* See k_exp.c for details.

*

* Input: ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7

* Output: 2**127 <= y < 2**128

*/

static float

__frexp_expf(float x, int *expt)

{

float exp_x;

uint32_t hx;

exp_x = expf(x - kln2);

GET_FLOAT_WORD(hx, exp_x);

*expt = (hx >> 23) - (0x7f + 127) + k;

SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23));

return (exp_x);

}

float

__ldexp_expf(float x, int expt)

{

float exp_x, scale;

int ex_expt;

exp_x = __frexp_expf(x, &ex_expt);

expt += ex_expt;

SET_FLOAT_WORD(scale, (0x7f + expt) << 23);

return (exp_x * scale);

}

float complex

__ldexp_cexpf(float complex z, int expt)

{

float x, y, exp_x, scale1, scale2;

int ex_expt, half_expt;

x = crealf(z);

y = cimagf(z);

exp_x = __frexp_expf(x, &ex_expt);

expt += ex_expt;

half_expt = expt / 2;

SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23);

half_expt = expt - half_expt;

SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);

return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,

sinf(y) * exp_x * scale1 * scale2));

}

#include <complex.h>

/*

* C99 specifies that complex numbers have the same representation as

* an array of two elements, where the first element is the real part

* and the second element is the imaginary part.

*/

typedef union {

float complex f;

float a[2];

} float_complex;

typedef union {

double complex f;

double a[2];

} double_complex;

typedef union {

long double complex f;

long double a[2];

} long_double_complex;

#define REALPART(z) ((z).a[0])

#define IMAGPART(z) ((z).a[1])

/*

* Inline functions that can be used to construct complex values.

*

* The C99 standard intends x+I*y to be used for this, but x+I*y is

* currently unusable in general since gcc introduces many overflow,

* underflow, sign and efficiency bugs by rewriting I*y as

* (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.

* In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted

* to -0.0+I*0.0.

*

* The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()

* to construct complex values. Compilers that conform to the C99

* standard require the following functions to avoid the above issues.

*/

#ifndef CMPLXF

static __inline float complex

CMPLXF(float x, float y)

{

double_complex z;

REALPART(z) = x;

IMAGPART(z) = y;

return (z.f);

}

#endif

#ifndef CMPLX

static __inline double complex

CMPLX(double x, double y)

{

double_complex z;

REALPART(z) = x;

IMAGPART(z) = y;

return (z.f);

}

#endif

#ifndef CMPLXL

static __inline long double complex

CMPLXL(long double x, long double y)

{

long_double_complex z;

REALPART(z) = x;

IMAGPART(z) = y;

return (z.f);

}

#endif

/* double precision kernel functions */

double __ldexp_exp(double,int);

#ifdef _COMPLEX_H_

double complex __ldexp_cexp(double complex,int);

#endif

/* float precision kernel functions */

float __ldexp_expf(float,int);

#ifdef _COMPLEX_H_

float complex __ldexp_cexpf(float complex,int);

#endif

/*-

* SPDX-License-Identifier: BSD-2-Clause-FreeBSD

* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>

* 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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 <sys/cdefs.h>

#include "math_private.h"

static const uint32_t

exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */

cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */

double x, y, exp_x;

uint32_t hx, hy, lx, ly;

x = creal(z);

y = cimag(z);

EXTRACT_WORDS(hy, ly, y);

hy &= 0x7fffffff;

/* cexp(x + I 0) = exp(x) + I 0 */

if ((hy | ly) == 0)

return (CMPLX(exp(x), y));

EXTRACT_WORDS(hx, lx, x);

/* cexp(0 + I y) = cos(y) + I sin(y) */

if (((hx & 0x7fffffff) | lx) == 0)

return (CMPLX(cos(y), sin(y)));

if (hy >= 0x7ff00000) {

if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {

/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */

return (CMPLX(y - y, y - y));

} else if (hx & 0x80000000) {

/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */

return (CMPLX(0.0, 0.0));

} else {

/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */

return (CMPLX(x, y - y));

}

}

if (hx >= exp_ovfl && hx <= cexp_ovfl) {

/*

* x is between 709.7 and 1454.3, so we must scale to avoid

* overflow in exp(x).

*/

return (__ldexp_cexp(z, 0));

} else {

/*

* Cases covered here:

* - x < exp_ovfl and exp(x) won't overflow (common case)

* - x > cexp_ovfl, so exp(x) * s overflows for all s > 0

* - x = +-Inf (generated by exp())

* - x = NaN (spurious inexact exception from y)

*/

exp_x = exp(x);

return (CMPLX(exp_x * cos(y), exp_x * sin(y)));

}

/*-

* SPDX-License-Identifier: BSD-2-Clause-FreeBSD

* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>

* 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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 <sys/cdefs.h>

#include "math_private.h"

static const uint32_t

exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */

cexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */

float x, y, exp_x;

uint32_t hx, hy;

x = crealf(z);

y = cimagf(z);

GET_FLOAT_WORD(hy, y);

hy &= 0x7fffffff;

/* cexp(x + I 0) = exp(x) + I 0 */

if (hy == 0)

return (CMPLXF(expf(x), y));

GET_FLOAT_WORD(hx, x);

/* cexp(0 + I y) = cos(y) + I sin(y) */

if ((hx & 0x7fffffff) == 0)

return (CMPLXF(cosf(y), sinf(y)));

if (hy >= 0x7f800000) {

if ((hx & 0x7fffffff) != 0x7f800000) {

/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */

return (CMPLXF(y - y, y - y));

} else if (hx & 0x80000000) {

/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */

return (CMPLXF(0.0, 0.0));

} else {

/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */

return (CMPLXF(x, y - y));

}

}

if (hx >= exp_ovfl && hx <= cexp_ovfl) {

/*

* x is between 88.7 and 192, so we must scale to avoid

* overflow in expf(x).

*/

return (__ldexp_cexpf(z, 0));

} else {

/*

* Cases covered here:

* - x < exp_ovfl and exp(x) won't overflow (common case)

* - x > cexp_ovfl, so exp(x) * s overflows for all s > 0

* - x = +-Inf (generated by exp())

* - x = NaN (spurious inexact exception from y)

*/

exp_x = expf(x);

return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));

}

DEF_STD(cexpf);

LDBL_MAYBE_UNUSED_CLONE(cexpf);

TESTS += cexp_test

/*-

* Copyright (c) 2008-2011 David Schultz <das@FreeBSD.org>

* 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, 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 AND CONTRIBUTORS ``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 OR CONTRIBUTORS 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.

*/

/*

* Tests for corner cases in cexp*().

*/

#include <sys/cdefs.h>

#include <sys/param.h>

#include <assert.h>

#include <complex.h>

#include <fenv.h>

#include <float.h>

#include <math.h>

#include <stdio.h>

#include "test-utils.h"

#pragma STDC FENV_ACCESS ON

#pragma STDC CX_LIMITED_RANGE OFF

/*

* Test that a function returns the correct value and sets the

* exception flags correctly. The exceptmask specifies which

* exceptions we should check. We need to be lenient for several

* reasons, but mainly because on some architectures it's impossible

* to raise FE_OVERFLOW without raising FE_INEXACT. In some cases,

* whether cexp() raises an invalid exception is unspecified.

*

* These are macros instead of functions so that assert provides more

* meaningful error messages.

*

* XXX The volatile here is to avoid gcc's bogus constant folding and work

* around the lack of support for the FENV_ACCESS pragma.

*/

#define test(func, z, result, exceptmask, excepts, checksign) do { \

volatile long double complex _d = z; \

assert(feclearexcept(FE_ALL_EXCEPT) == 0); \

assert(cfpequal_cs((func)(_d), (result), (checksign))); \

assert(((void)(func), fetestexcept(exceptmask) == (excepts))); \

} while (0)

#define test_c(t, func, z, result, exceptmask, excepts, checksign) do { \

volatile t complex r = result; \

test(func, z, r, exceptmask, excepts, checksign); \

} while (0)

#define test_f(func, z, result, exceptmask, excepts, checksign) do { \

test_c(float, func, z, result, exceptmask, excepts, checksign); \

} while (0)

/* Test within a given tolerance. */

#define test_tol(func, z, result, tol) do { \

volatile long double complex _d = z; \

assert(cfpequal_tol((func)(_d), (result), (tol), \

FPE_ABS_ZERO | CS_BOTH)); \

} while (0)

/* Test all the functions that compute cexp(x). */

#define testall(x, result, exceptmask, excepts, checksign) do { \

test(cexp, x, result, exceptmask, excepts, checksign); \

test_f(cexpf, x, result, exceptmask, excepts, checksign); \

} while (0)

/*

* Test all the functions that compute cexp(x), within a given tolerance.

* The tolerance is specified in ulps.

*/

#define testall_tol(x, result, tol) do { \

test_tol(cexp, x, result, tol * DBL_ULP()); \

test_tol(cexpf, x, result, tol * FLT_ULP()); \

} while (0)

/* Various finite non-zero numbers to test. */

static const float finites[] =

{ -42.0e20, -1.0, -1.0e-10, -0.0, 0.0, 1.0e-10, 1.0, 42.0e20 };

/* Tests for 0 */

static void

test_zero(void)

{

/* cexp(0) = 1, no exceptions raised */

testall(0.0, 1.0, ALL_STD_EXCEPT, 0, 1);

testall(-0.0, 1.0, ALL_STD_EXCEPT, 0, 1);

testall(CMPLXL(0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);

testall(CMPLXL(-0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1);

}

/*

* Tests for NaN. The signs of the results are indeterminate unless the

* imaginary part is 0.

*/

static void

test_nan(void)

{

unsigned i;

/* cexp(x + NaNi) = NaN + NaNi and optionally raises invalid */

/* cexp(NaN + yi) = NaN + NaNi and optionally raises invalid (|y|>0) */

for (i = 0; i < nitems(finites); i++) {

printf("# Run %d..\n", i);

testall(CMPLXL(finites[i], NAN), CMPLXL(NAN, NAN),

ALL_STD_EXCEPT & ~FE_INVALID, 0, 0);

if (finites[i] == 0.0)

continue;

/* XXX FE_INEXACT shouldn't be raised here */

testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN),

ALL_STD_EXCEPT & ~(FE_INVALID | FE_INEXACT), 0, 0);

}

/* cexp(NaN +- 0i) = NaN +- 0i */

testall(CMPLXL(NAN, 0.0), CMPLXL(NAN, 0.0), ALL_STD_EXCEPT, 0, 1);

testall(CMPLXL(NAN, -0.0), CMPLXL(NAN, -0.0), ALL_STD_EXCEPT, 0, 1);

/* cexp(inf + NaN i) = inf + nan i */

testall(CMPLXL(INFINITY, NAN), CMPLXL(INFINITY, NAN),

ALL_STD_EXCEPT, 0, 0);

/* cexp(-inf + NaN i) = 0 */

testall(CMPLXL(-INFINITY, NAN), CMPLXL(0.0, 0.0),

ALL_STD_EXCEPT, 0, 0);

/* cexp(NaN + NaN i) = NaN + NaN i */

testall(CMPLXL(NAN, NAN), CMPLXL(NAN, NAN),

ALL_STD_EXCEPT, 0, 0);

}

static void

test_inf(void)

{

unsigned i;

/* cexp(x + inf i) = NaN + NaNi and raises invalid */

for (i = 0; i < nitems(finites); i++) {

printf("# Run %d..\n", i);

testall(CMPLXL(finites[i], INFINITY), CMPLXL(NAN, NAN),

ALL_STD_EXCEPT, FE_INVALID, 1);

}

/* cexp(-inf + yi) = 0 * (cos(y) + sin(y)i) */

/* XXX shouldn't raise an inexact exception */

testall(CMPLXL(-INFINITY, M_PI_4), CMPLXL(0.0, 0.0),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(-INFINITY, 3 * M_PI_4), CMPLXL(-0.0, 0.0),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(-INFINITY, 5 * M_PI_4), CMPLXL(-0.0, -0.0),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(-INFINITY, 7 * M_PI_4), CMPLXL(0.0, -0.0),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(-INFINITY, 0.0), CMPLXL(0.0, 0.0),

ALL_STD_EXCEPT, 0, 1);

testall(CMPLXL(-INFINITY, -0.0), CMPLXL(0.0, -0.0),

ALL_STD_EXCEPT, 0, 1);

/* cexp(inf + yi) = inf * (cos(y) + sin(y)i) (except y=0) */

/* XXX shouldn't raise an inexact exception */

testall(CMPLXL(INFINITY, M_PI_4), CMPLXL(INFINITY, INFINITY),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(INFINITY, 3 * M_PI_4), CMPLXL(-INFINITY, INFINITY),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(INFINITY, 5 * M_PI_4), CMPLXL(-INFINITY, -INFINITY),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

testall(CMPLXL(INFINITY, 7 * M_PI_4), CMPLXL(INFINITY, -INFINITY),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

/* cexp(inf + 0i) = inf + 0i */

testall(CMPLXL(INFINITY, 0.0), CMPLXL(INFINITY, 0.0),

ALL_STD_EXCEPT, 0, 1);

testall(CMPLXL(INFINITY, -0.0), CMPLXL(INFINITY, -0.0),

ALL_STD_EXCEPT, 0, 1);

}

static void

test_reals(void)

{

unsigned i;

for (i = 0; i < nitems(finites); i++) {

/* XXX could check exceptions more meticulously */

printf("# Run %d..\n", i);

test(cexp, CMPLXL(finites[i], 0.0),

CMPLXL(exp(finites[i]), 0.0),

FE_INVALID | FE_DIVBYZERO, 0, 1);

test(cexp, CMPLXL(finites[i], -0.0),

CMPLXL(exp(finites[i]), -0.0),

FE_INVALID | FE_DIVBYZERO, 0, 1);

test_f(cexpf, CMPLXL(finites[i], 0.0),

CMPLXL(expf(finites[i]), 0.0),

FE_INVALID | FE_DIVBYZERO, 0, 1);

test_f(cexpf, CMPLXL(finites[i], -0.0),

CMPLXL(expf(finites[i]), -0.0),

FE_INVALID | FE_DIVBYZERO, 0, 1);

}

}

static void

test_imaginaries(void)

{

unsigned i;

for (i = 0; i < nitems(finites); i++) {

printf("# Run %d..\n", i);

test(cexp, CMPLXL(0.0, finites[i]),

CMPLXL(cos(finites[i]), sin(finites[i])),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

test(cexp, CMPLXL(-0.0, finites[i]),

CMPLXL(cos(finites[i]), sin(finites[i])),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

test_f(cexpf, CMPLXL(0.0, finites[i]),

CMPLXL(cosf(finites[i]), sinf(finites[i])),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

test_f(cexpf, CMPLXL(-0.0, finites[i]),

CMPLXL(cosf(finites[i]), sinf(finites[i])),

ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1);

}

}

static void

test_small(void)

{

static const double tests[] = {

/* csqrt(a + bI) = x + yI */

/* a b x y */

1.0, M_PI_4, M_SQRT2 * 0.5 * M_E, M_SQRT2 * 0.5 * M_E,

-1.0, M_PI_4, M_SQRT2 * 0.5 / M_E, M_SQRT2 * 0.5 / M_E,

2.0, M_PI_2, 0.0, M_E * M_E,

M_LN2, M_PI, -2.0, 0.0,

};

double a, b;

double x, y;

unsigned i;

for (i = 0; i < nitems(tests); i += 4) {

printf("# Run %d..\n", i);

a = tests[i];

b = tests[i + 1];

x = tests[i + 2];

y = tests[i + 3];

test_tol(cexp, CMPLXL(a, b), CMPLXL(x, y), 3 * DBL_ULP());

/* float doesn't have enough precision to pass these tests */

if (x == 0 || y == 0)

continue;

test_tol(cexpf, CMPLXL(a, b), CMPLXL(x, y), 1 * FLT_ULP());

}

}

/* Test inputs with a real part r that would overflow exp(r). */

static void

test_large(void)

{

test_tol(cexp, CMPLXL(709.79, 0x1p-1074),

CMPLXL(INFINITY, 8.94674309915433533273e-16), DBL_ULP());

test_tol(cexp, CMPLXL(1000, 0x1p-1074),

CMPLXL(INFINITY, 9.73344457300016401328e+110), DBL_ULP());

test_tol(cexp, CMPLXL(1400, 0x1p-1074),

CMPLXL(INFINITY, 5.08228858149196559681e+284), DBL_ULP());

test_tol(cexp, CMPLXL(900, 0x1.23456789abcdep-1020),

CMPLXL(INFINITY, 7.42156649354218408074e+83), DBL_ULP());

test_tol(cexp, CMPLXL(1300, 0x1.23456789abcdep-1020),

CMPLXL(INFINITY, 3.87514844965996756704e+257), DBL_ULP());

test_tol(cexpf, CMPLXL(88.73, 0x1p-149),

CMPLXL(INFINITY, 4.80265603e-07), 2 * FLT_ULP());

test_tol(cexpf, CMPLXL(90, 0x1p-149),

CMPLXL(INFINITY, 1.7101492622e-06f), 2 * FLT_ULP());

test_tol(cexpf, CMPLXL(192, 0x1p-149),

CMPLXL(INFINITY, 3.396809344e+38f), 2 * FLT_ULP());

test_tol(cexpf, CMPLXL(120, 0x1.234568p-120),

CMPLXL(INFINITY, 1.1163382522e+16f), 2 * FLT_ULP());

test_tol(cexpf, CMPLXL(170, 0x1.234568p-120),

CMPLXL(INFINITY, 5.7878851079e+37f), 2 * FLT_ULP());

}

int

main(void)

{

printf("1..7\n");

test_zero();

printf("ok 1 - cexp zero\n");

test_nan();

printf("ok 2 - cexp nan\n");

test_inf();

printf("ok 3 - cexp inf\n");

test_reals();

printf("ok 4 - cexp reals\n");

test_imaginaries();

printf("ok 5 - cexp imaginaries\n");

test_small();

printf("ok 6 - cexp small\n");

test_large();

printf("ok 7 - cexp large\n");

return (0);

}