| s_ctanh.c (226600) | s_ctanh.c (275819) |
|---|---|
| 1/*- 2 * Copyright (c) 2011 David Schultz 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 --- 50 unchanged lines hidden (view full) --- 59 * Modifications: 60 * 61 * I omitted the original algorithm's handling of overflow in tan(x) after 62 * verifying with nearpi.c that this can't happen in IEEE single or double 63 * precision. I also handle large x differently. 64 */ 65 66#include <sys/cdefs.h> | 1/*- 2 * Copyright (c) 2011 David Schultz 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 --- 50 unchanged lines hidden (view full) --- 59 * Modifications: 60 * 61 * I omitted the original algorithm's handling of overflow in tan(x) after 62 * verifying with nearpi.c that this can't happen in IEEE single or double 63 * precision. I also handle large x differently. 64 */ 65 66#include <sys/cdefs.h> |
| 67__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanh.c 226600 2011-10-21 06:30:16Z das $"); | 67__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanh.c 275819 2014-12-16 09:21:56Z ed $"); |
| 68 69#include <complex.h> 70#include <math.h> 71 72#include "math_private.h" 73 74double complex 75ctanh(double complex z) --- 21 unchanged lines hidden (view full) --- 97 * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite 98 * 99 * The imaginary part of the sign is unspecified. This special 100 * case is only needed to avoid a spurious invalid exception when 101 * y is infinite. 102 */ 103 if (ix >= 0x7ff00000) { 104 if ((ix & 0xfffff) | lx) /* x is NaN */ | 68 69#include <complex.h> 70#include <math.h> 71 72#include "math_private.h" 73 74double complex 75ctanh(double complex z) --- 21 unchanged lines hidden (view full) --- 97 * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite 98 * 99 * The imaginary part of the sign is unspecified. This special 100 * case is only needed to avoid a spurious invalid exception when 101 * y is infinite. 102 */ 103 if (ix >= 0x7ff00000) { 104 if ((ix & 0xfffff) | lx) /* x is NaN */ |
| 105 return (cpack(x, (y == 0 ? y : x * y))); | 105 return (CMPLX(x, (y == 0 ? y : x * y))); |
| 106 SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ | 106 SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ |
| 107 return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); | 107 return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); |
| 108 } 109 110 /* 111 * ctanh(x + i NAN) = NaN + i NaN 112 * ctanh(x +- i Inf) = NaN + i NaN 113 */ 114 if (!isfinite(y)) | 108 } 109 110 /* 111 * ctanh(x + i NAN) = NaN + i NaN 112 * ctanh(x +- i Inf) = NaN + i NaN 113 */ 114 if (!isfinite(y)) |
| 115 return (cpack(y - y, y - y)); | 115 return (CMPLX(y - y, y - y)); |
| 116 117 /* 118 * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the 119 * approximation sinh^2(huge) ~= exp(2*huge) / 4. 120 * We use a modified formula to avoid spurious overflow. 121 */ 122 if (ix >= 0x40360000) { /* x >= 22 */ 123 double exp_mx = exp(-fabs(x)); | 116 117 /* 118 * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the 119 * approximation sinh^2(huge) ~= exp(2*huge) / 4. 120 * We use a modified formula to avoid spurious overflow. 121 */ 122 if (ix >= 0x40360000) { /* x >= 22 */ 123 double exp_mx = exp(-fabs(x)); |
| 124 return (cpack(copysign(1, x), | 124 return (CMPLX(copysign(1, x), |
| 125 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 126 } 127 128 /* Kahan's algorithm */ 129 t = tan(y); 130 beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 131 s = sinh(x); 132 rho = sqrt(1 + s * s); /* = cosh(x) */ 133 denom = 1 + beta * s * s; | 125 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 126 } 127 128 /* Kahan's algorithm */ 129 t = tan(y); 130 beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 131 s = sinh(x); 132 rho = sqrt(1 + s * s); /* = cosh(x) */ 133 denom = 1 + beta * s * s; |
| 134 return (cpack((beta * rho * s) / denom, t / denom)); | 134 return (CMPLX((beta * rho * s) / denom, t / denom)); |
| 135} 136 137double complex 138ctan(double complex z) 139{ 140 141 /* ctan(z) = -I * ctanh(I * z) */ | 135} 136 137double complex 138ctan(double complex z) 139{ 140 141 /* ctan(z) = -I * ctanh(I * z) */ |
| 142 z = ctanh(cpack(-cimag(z), creal(z))); 143 return (cpack(cimag(z), -creal(z))); | 142 z = ctanh(CMPLX(-cimag(z), creal(z))); 143 return (CMPLX(cimag(z), -creal(z))); |
| 144} | 144} |