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