1 1.1 drochner /* $NetBSD: cephes_subrf.c,v 1.1 2007/08/20 16:01:34 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_subrf.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 _cchshf(float x, float *c, float *s) 41 1.1 drochner { 42 1.1 drochner float e, ei; 43 1.1 drochner 44 1.1 drochner if (fabsf(x) <= 0.5f) { 45 1.1 drochner *c = coshf(x); 46 1.1 drochner *s = sinhf(x); 47 1.1 drochner } else { 48 1.1 drochner e = expf(x); 49 1.1 drochner ei = 0.5f / e; 50 1.1 drochner e = 0.5f * 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.140625; 60 1.1 drochner static const double DP2 = 9.67502593994140625E-4; 61 1.1 drochner static const double DP3 = 1.509957990978376432E-7; 62 1.1 drochner #define MACHEPF 3.0e-8 63 1.1 drochner 64 1.1 drochner float 65 1.1 drochner _redupif(float x) 66 1.1 drochner { 67 1.1 drochner float t; 68 1.1 drochner long i; 69 1.1 drochner 70 1.1 drochner t = x / (float)M_PI; 71 1.1 drochner if (t >= 0.0f) 72 1.1 drochner t += 0.5f; 73 1.1 drochner else 74 1.1 drochner t -= 0.5f; 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 float 85 1.1 drochner _ctansf(float complex z) 86 1.1 drochner { 87 1.1 drochner float f, x, x2, y, y2, rn, t, d; 88 1.1 drochner 89 1.1 drochner x = fabsf(2.0f * crealf(z)); 90 1.1 drochner y = fabsf(2.0f * cimagf(z)); 91 1.1 drochner 92 1.1 drochner x = _redupif(x); 93 1.1 drochner 94 1.1 drochner x = x * x; 95 1.1 drochner y = y * y; 96 1.1 drochner x2 = 1.0f; 97 1.1 drochner y2 = 1.0f; 98 1.1 drochner f = 1.0f; 99 1.1 drochner rn = 0.0f; 100 1.1 drochner d = 0.0f; 101 1.1 drochner do { 102 1.1 drochner rn += 1.0f; 103 1.1 drochner f *= rn; 104 1.1 drochner rn += 1.0f; 105 1.1 drochner f *= rn; 106 1.1 drochner x2 *= x; 107 1.1 drochner y2 *= y; 108 1.1 drochner t = y2 + x2; 109 1.1 drochner t /= f; 110 1.1 drochner d += t; 111 1.1 drochner 112 1.1 drochner rn += 1.0f; 113 1.1 drochner f *= rn; 114 1.1 drochner rn += 1.0f; 115 1.1 drochner f *= rn; 116 1.1 drochner x2 *= x; 117 1.1 drochner y2 *= y; 118 1.1 drochner t = y2 - x2; 119 1.1 drochner t /= f; 120 1.1 drochner d += t; 121 1.1 drochner } while (fabsf(t/d) > MACHEPF); 122 1.1 drochner return d; 123 1.1 drochner } 124