1/* $NetBSD$ */
2
3/*-
4 * Copyright (c) 2007 The NetBSD Foundation, Inc.
5 * All rights reserved.
6 *
7 * This code is derived from software written by Stephen L. Moshier.
8 * It is redistributed by the NetBSD Foundation by permission of the author.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 *    notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 *    notice, this list of conditions and the following disclaimer in the
17 *    documentation and/or other materials provided with the distribution.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
20 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
21 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
23 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 */
31
32#include "../src/namespace.h"
33#include <complex.h>
34#include <math.h>
35#include "cephes_subr.h"
36
37/* calculate cosh and sinh */
38
39void
40_cchsh(double x, double *c, double *s)
41{
42	double e, ei;
43
44	if (fabs(x) <= 0.5) {
45		*c = cosh(x);
46		*s = sinh(x);
47	} else {
48		e = exp(x);
49		ei = 0.5 / e;
50		e = 0.5 * e;
51		*s = e - ei;
52		*c = e + ei;
53	}
54}
55
56/* Program to subtract nearest integer multiple of PI */
57
58/* extended precision value of PI: */
59static const double DP1 = 3.14159265160560607910E0;
60static const double DP2 = 1.98418714791870343106E-9;
61static const double DP3 = 1.14423774522196636802E-17;
62#define MACHEP 1.1e-16
63
64double
65_redupi(double x)
66{
67	double t;
68	long i;
69
70	t = x / M_PI;
71	if (t >= 0.0)
72		t += 0.5;
73	else
74		t -= 0.5;
75
76	i = t;	/* the multiple */
77	t = i;
78	t = ((x - t * DP1) - t * DP2) - t * DP3;
79	return t;
80}
81
82/* Taylor series expansion for cosh(2y) - cos(2x) */
83
84double
85_ctans(double complex z)
86{
87	double f, x, x2, y, y2, rn, t;
88	double d;
89
90	x = fabs(2.0 * creal(z));
91	y = fabs(2.0 * cimag(z));
92
93	x = _redupi(x);
94
95	x = x * x;
96	y = y * y;
97	x2 = 1.0;
98	y2 = 1.0;
99	f = 1.0;
100	rn = 0.0;
101	d = 0.0;
102	do {
103		rn += 1.0;
104		f *= rn;
105		rn += 1.0;
106		f *= rn;
107		x2 *= x;
108		y2 *= y;
109		t = y2 + x2;
110		t /= f;
111		d += t;
112
113		rn += 1.0;
114		f *= rn;
115		rn += 1.0;
116		f *= rn;
117		x2 *= x;
118		y2 *= y;
119		t = y2 - x2;
120		t /= f;
121		d += t;
122	} while (fabs(t/d) > MACHEP);
123	return d;
124}
125