/* * Copyright (c) 2002 Apple Computer, Inc. All rights reserved. * * @APPLE_LICENSE_HEADER_START@ * * The contents of this file constitute Original Code as defined in and * are subject to the Apple Public Source License Version 1.1 (the * "License"). You may not use this file except in compliance with the * License. Please obtain a copy of the License at * http://www.apple.com/publicsource and read it before using this file. * * This Original Code and all software distributed under the License are * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES, * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the * License for the specific language governing rights and limitations * under the License. * * @APPLE_LICENSE_HEADER_END@ */ /**************************************************************************** ** File: complex.c ** ** Contains: C source code for implementations of floating-point ** (double) complex functions defined in header file ** "complex.h" for PowerPC Macintoshes in native mode. ** Transcendental function algorithms are based on the ** paper "Branch Cuts for Complex Elementary Functions" ** by W. Kahan, May 17, 1987, and on Pascal and C source ** code for the SANE 80-/96-bit extended type by Kenton ** Hanson and Paul Finlayson, respectively. ** ** ** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838 ** ** Copyright: c 1987-1993 by Apple Computer, Inc., all rights reserved. ** ** Change History (most recent first): ** ** 25 Aug 93 ali Changed clog to cLog to avoid clashing with the ** stream i/o definition clog. ** 14 Jul 93 ali Added #pragma fenv_access on ** 22 Feb 93 ali Added a nomaf #pragma. ** 05 Feb 93 JPO Modified calls to feclearexcept and feraiseexcept ** to reflect changes in "fenv.h". ** 18 Dec 92 JPO First created. ** ****************************************************************************/ #pragma option nomaf #pragma fenv_access on #include "math.h" #include "complex.h" #include "fenv.h" #define Real(z) (__real__ z) #define Imag(z) (__imag__ z) /**************************************************************************** CONSTANTS used by complex functions ****************************************************************************/ #if defined(__BIG_ENDIAN__) #define HEXDOUBLE(hi, lo) { { hi, lo } } #elif defined(__LITTLE_ENDIAN__) #define HEXDOUBLE(hi, lo) { { lo, hi } } #else #error Unknown endianness #endif static const union { /* +infinity */ long int ival[2]; double dval; } FPKINF = HEXDOUBLE(0x7ff00000,0x00000000); static const union { /* sqrt(2.0) */ long int ival[2]; double dval; } FPKSQT2 = HEXDOUBLE(0x3ff6a09e,0x667f3bcd); static const union { /* sqrt(2.0) + 1.0 to double */ long int ival[2]; double dval; } FPKR2P1 = HEXDOUBLE(0x4003504f,0x333f9de6); static const union { /* sqrt(2.0) + 1.0 - FPKR2P1 to double */ long int ival[2]; double dval; } FPKT2P1 = HEXDOUBLE(0x3ca21165,0xf626cdd6); static const union { /* pi/2.0 */ long int ival[2]; double dval; } FPKPI2 = HEXDOUBLE(0x3ff921fb,0x54442d18); static const union { /* pi */ long int ival[2]; double dval; } FPKPI = HEXDOUBLE(0x400921fb,0x54442d18); static const union { /* underflow threshold / round threshold */ long int ival[2]; double dval; } FPKLOVERE = HEXDOUBLE(0x03600000,0x00000000); static const union { /* sqrt(0.5) */ long int ival[2]; double dval; } FPKSQRTHALF = HEXDOUBLE(0x3fe6a09e,0x667f3bcd); static const union { /* log(2.0) */ long int ival[2]; double dval; } FPKLOG2 = HEXDOUBLE(0x3fe62e42,0xfefa39ef); static const union { /* exp(709.0) */ long int ival[2]; double dval; } FPKEXP709 = HEXDOUBLE(0x7fdd422d,0x2be5dc9b); static const union { /* asinh(nextafterd(+infinity,0.0))/4.0 */ long int ival[2]; double dval; } FPKASINHOM4 = HEXDOUBLE(0x406633ce,0x8fb9f87e); static const union { /* 2.0^(-33) */ long int ival[2]; double dval; } FPK2M33 = HEXDOUBLE(0x3de00000,0x00000000); static const union { /* 2.0^54 */ long int ival[2]; double dval; } FPK2P54 = HEXDOUBLE(0x43500000,0x00000000); #if 0 static const union { /* 2.0^-(1022) */ long int ival[2]; double dval; } FPKTINY = HEXDOUBLE(0x00100000,0x00000000); #endif static const union { /* sqrt(nextafterd(+infinity,0.0))/4.0 */ long int ival[2]; double dval; } FPKTHETA = HEXDOUBLE(0x5fcfffff,0xffffffff); static const union { /* 4.0/sqrt(nextafterd(+infinity,0.0)) */ long int ival[2]; double dval; } FPKRHO = HEXDOUBLE(0x20100000,0x00000000); /**************************************************************************** double complex xdivc(double x, double complex y) returns (real x) / (complex y). CONSTANTS: FPKINF = +infinity ****************************************************************************/ static double complex xdivc( double x, double complex y ) /* returns (real x) / (complex y) */ { double complex z; double r, denom; if ( fabs(Real(y)) >= fabs(Imag(y)) ) { /* |Real(y)| >= |Imag(y)| */ if (fabs(Real(y)) == FPKINF.dval) { /* Imag(y) and Real(y) are infinite */ Real(z) = copysign(0.0,Real(y)); Imag(z) = copysign(0.0,-Imag(y)); } else { /* |Real(y)| >= finite |Imag(y)| */ r = Imag(y)/Real(y); denom = Real(y) + Imag(y)*r; Real(z) = x/denom; Imag(z) = (-x*r)/denom; } } else { /* |Real(y)| !>= |Imag(y)| */ r = Real(y)/Imag(y); denom = r*Real(y) + Imag(y); Real(z) = (r*x)/denom; Imag(z) = -x/denom; } return z; } #ifdef __ppc__ static double __cabs ( double complex z ) { double a,b,s,t; fenv_t env; int clre,clim,ifoo; clre = fpclassify(Real(z)); clim = fpclassify(Imag(z)); if ((clre < FP_NORMAL) || (clim < FP_NORMAL)) { return (fabs(Real(z)) + fabs(Imag(z))); /* Real(z) or Imag(z) is NaN, INF, or zero */ } else { /* both components of z are finite, nonzero */ ifoo = feholdexcept(&env); /* save environment, clear flags */ ifoo = fesetround(FE_TONEAREST); /* set default rounding */ a = fabs(Real(z)); /* work with absolute values */ b = fabs(Imag(z)); s = 0.0; if (a < b) { /* order a >= b */ t = a; a = b; b = t; } t = a - b; /* magnitude difference */ if (t != a) { /* b not negligible relative to a */ if (t > b) { /* a - b > b */ s = a/b; s += sqrt(1.0 + s*s); } else { /* a - b <= b */ s = t/b; t = (2.0 + s)*s; s = ((FPKT2P1.dval+t/(FPKSQT2.dval+sqrt(2.0+t)))+s)+FPKR2P1.dval; } s = b/s; /* may spuriously underflow */ feclearexcept(FE_UNDERFLOW); } feupdateenv(&env); /* restore environment */ return (a + s); /* deserved overflow occurs here */ } /* finite, nonzero case */ } // XXX Revert __cabs call to cabs and discard __cabs above post-10.3 (chicken/egg situation) float cabsf( float complex z ){ return (float) __cabs((double complex) z); } #endif /**************************************************************************** double carg(double complex z) returns the argument (in radians) of its complex argument z. The algorithm is from Kahan's paper. The argument of a complex number z = x + i*y is defined as atan2(y,x) for finite x and y. CONSTANTS: FPKPI2 = pi/2.0 to double precision FPKPI = pi to double precision Calls: fpclassify, copysign, fabs, atan ****************************************************************************/ double carg ( double complex z ) { double a,b,argr; int clre,clim; a = Real(z); b = Imag(z); clre = fpclassify(a); clim = fpclassify(b); if ((clre == FP_ZERO) && (clim == FP_ZERO)) { /* zero real and imag parts */ a = copysign(1.0,a); } if ((clre == FP_INFINITE) && (clim == FP_INFINITE)) { /* both parts INF */ a = copysign(1.0,a); b = copysign(1.0,b); } if (fabs(b) > fabs(a)) /* |imag| > |real| */ argr = copysign(FPKPI2.dval,b) - atan(a/b); else { if (a < 0.0) /* |real| >= |imag| */ argr = copysign(FPKPI.dval,b) + atan(b/a); else argr = atan(b/a); } return argr; } /**************************************************************************** static double cssqs(double complex z, long int *k) returns rho = |z/(2^*k)|^2 with scale factor *k set to avoid overflow/underflow. Algorithm is from the Kahan paper. CONSTANTS: FPKINF = +infinity FPKLOVERE = (double underflow threshold)/(double round threshold) = (4.0*nextafterd(1.0,0.0)/nextafterd(INF,0.0)) /(1.0 - nextafterd(1.0,0) Calls: fabs, logb, scalb, fmax, feholdexcept, fetestexcept, feclearexcept, and feupdateenv. Called by: csqrt and clog. ****************************************************************************/ static double cssqs ( double complex z, long int *k) { double a,b,rho; fenv_t env; long int iscale; int ifoo; iscale = 0; a = fabs(Real(z)); b = fabs(Imag(z)); ifoo = feholdexcept(&env); /* save environment, clr flags */ rho = a*a + b*b; /* preliminary result */ if ((a == FPKINF.dval) || (b == FPKINF.dval)) { rho = FPKINF.dval; /* +INF result if Real(z) or Imag(z) is infinite */ } else if (fetestexcept(FE_OVERFLOW) || (fetestexcept(FE_UNDERFLOW) && (rho < FPKLOVERE.dval))) { iscale = (long int) logb(fmax(a,b)); /* scaling necessary */ a = scalb(a,-iscale); b = scalb(b,-iscale); rho = a*a + b*b; /* re-calculate scaled square magnitude */ } feclearexcept(FE_OVERFLOW + FE_UNDERFLOW); feupdateenv(&env); /* restore environment */ *k = iscale; /* store scale value */ return (rho); } /**************************************************************************** double complex csqrt(double complex z) returns the complex square root of its argument. The algorithm, which is from the Kahan paper, uses the following identities: sqrt(x + i*y) = sqrt((|z| + Real(z))/2) + i*sqrt((|z| - Real(z))/2) and sqrt(x - i*y) = sqrt((|z| + Real(z))/2) - i*sqrt((|z| - Real(z))/2), where y is positive and x may be positive or negative. CONSTANTS: FPKINF = +infinity Calls: cssqs, scalb, fabs, sqrt, copysign. ****************************************************************************/ double complex csqrt ( double complex z ) { double rho,x,y; double complex w; long int k; rho = cssqs(z,&k); /* scaled square magnitude */ if (Real(z) == Real(z)) rho = scalb(fabs(Real(z)),-k) + sqrt(rho); /* scaled |Real(z)| + |z| */ if (k%2) /* k is odd */ k = (k-1)/2; else { /* k is even */ k = k/2 - 1; rho = rho + rho; } rho = scalb(sqrt(rho),k); /* sqrt((|Real(z)| + |z|)/2) */ x = rho; y = Imag(z); if (rho != 0) { if (fabs(y) != FPKINF.dval) y = (y/rho)*0.5; /* signbit(Imag(z))*sqrt((|z|-|Real(z)|)/2 */ if (Real(z) < 0) { x = fabs(y); y = copysign(rho,Imag(z)); } } Real(w) = x; Imag(w) = y; return w; } /**************************************************************************** double complex clog(double complex z) returns the complex natural logarithm of its argument. The algorithm, which is from the Kahan paper, avoids spurious underflow or overflow. CONSTANTS: FPKSQRTHALF = sqrt(0.5) to double precision FPKLOG2 = log(2.0) to double precision Calls: fabs, cssqs, log1p, log, carg. ****************************************************************************/ double complex clog ( double complex z ) { double rho,dmax,dmin,temp; double complex w; long int k; dmax = fabs(Real(z)); /* order real and imaginary parts of z by magnitude */ dmin = fabs(Imag(z)); if (dmax < dmin) { temp = dmax; dmax = dmin; dmin = temp; } rho = cssqs(z,&k); /* scaled |z*z| */ if ((k == 0) && (dmax > FPKSQRTHALF.dval) && ((dmax <= 1.25) || (rho < 3.0))) Real(w) = log1p((dmax - 1.0)*(dmax + 1.0) + dmin*dmin)*0.5; /* |z| near 1.0 */ else Real(w) = log(rho)*0.5 + k*FPKLOG2.dval; /* more naive approximation */ Imag(w) = carg(z); /* imaginary part of logarithm */ return w; } /**************************************************************************** static double coshmul(double x, double y) returns y*cosh(x) while avoiding spurious overflow and invalid exceptions. CONSTANTS: FPKEXP709 = exp(709.0) in double precision Calls: cosh, exp, fpclassify, fabs, and scalb. Called by: csin, ccos, csinh, and ccosh. ****************************************************************************/ static double coshmul ( double x, double y ) { double absx, result; absx = fabs(x); if (absx <= 709.0) { /* cosh(x) is finite */ return y*cosh(x); } else if (fpclassify(x) < FP_ZERO) { /* x is NaN or infinite */ return (y*absx); } else if (absx > 1460.0) { /* probable overflow case */ return (scalb(y,2100)); } else { /* cosh(x) overflows but y*cosh(x) may not */ result = (0.5 * FPKEXP709.dval) * y; /* initialize result to cosh(709) */ absx -= 709.0; /* reduce exponential argument */ while (absx > 709.0) { /* exponential reduction loop */ result *= FPKEXP709.dval; absx -= 709.0; } return (result*exp(absx)); /* final multiplication */ } } /**************************************************************************** static double sinhmul(double x, double y) returns y*sinh(x) while avoiding spurious overflow and invalid exceptions. CONSTANTS: FPKEXP709 = exp(709.0) in double precision Calls: sinh, exp, fpclassify, fabs, and scalb. Called by: csin, ccos, csinh, and ccosh. ****************************************************************************/ static double sinhmul ( double x, double y ) { double absx, result; absx = fabs(x); if (absx <= 709.0) { /* sinh(x) is finite */ return y*sinh(x); } else if (fpclassify(x) < FP_ZERO) { /* x is NaN or infinite */ return (y*x); } else if (absx > 1460.0) { /* probable overflow case */ return (scalb(y,2100)); } else { /* sinh(x) overflows but y*sinh(x) may not */ result = (0.5*FPKEXP709.dval)*y; /* initialize result to |sinh(709)| */ absx -= 709.0; /* reduce exponential argument */ if (signbit(x) != 0) result = -result; /* take care of sign of result */ while (absx > 709.0) { /* exponential reduction loop */ result *= FPKEXP709.dval; absx -= 709.0; } return (result*exp(absx)); /* final multiplication */ } } /**************************************************************************** double complex csin(double complex z) returns the complex sine of its argument. The algorithm is based upon the identity: sin(x + i*y) = sin(x)*cosh(y) + i*cos(x)*sinh(y). Signaling of spurious overflows, underflows, and invalids is avoided where possible. Calls: sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex csin ( double complex z ) { fenv_t env; int ifoo; double sinval, cosval; double complex w; ifoo = feholdexcept(&env); /* save environment, clear flags */ sinval = sin(Real(z)); /* sine of real part of argument */ cosval = cos(Real(z)); /* cosine of real part of argument */ Real(w) = coshmul(Imag(z),sinval); /* real result = sinval*cosh(Imag(z)) */ Imag(w) = sinhmul(Imag(z),cosval); /* imag result = cosval*sinh(Imag(z)) */ feclearexcept(FE_UNDERFLOW); /* underflows don't matter */ feupdateenv(&env); /* update environment */ return w; } /**************************************************************************** double complex ccos(double complex z) returns the complex cosine of its argument. The algorithm is based upon the identity: cos(x + i*y) = cos(x)*cosh(y) - i*sin(x)*sinh(y). Signaling of spurious overflows, underflows, and invalids is avoided where possible. Calls: sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex ccos ( double complex z ) { fenv_t env; int ifoo; double sinval, cosval; double complex w; ifoo = feholdexcept(&env); /* save environment, clear flags */ sinval = sin(Real(z)); /* sine of real part of argument */ cosval = cos(Real(z)); /* cosine of real part of argument */ Real(w) = coshmul(Imag(z),cosval); /* real result = cosval*cosh(Imag(z)) */ Imag(w) = sinhmul(Imag(z),-sinval); /* imag result = -sinval*sinh(Imag(z)) */ feclearexcept(FE_UNDERFLOW); /* underflows don't matter */ feupdateenv(&env); /* update environment */ return w; } /**************************************************************************** double complex csinh(double complex z) returns the complex hyperbolic sine of its argument. The algorithm is based upon the identity: sinh(x + i*y) = cos(y)*sinh(x) + i*sin(y)*cosh(x). Signaling of spurious overflows, underflows, and invalids is avoided where possible. Calls: sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex csinh ( double complex z ) { fenv_t env; int ifoo; double sinval, cosval; double complex w; ifoo = feholdexcept(&env); /* save environment, clear flags */ sinval = sin(Imag(z)); /* sine of imaginary part of argument */ cosval = cos(Imag(z)); /* cosine of imaginary part of argument */ Real(w) = sinhmul(Real(z),cosval); /* real result = cosval*sinh(Real(z)) */ Imag(w) = coshmul(Real(z),sinval); /* imag result = sinval*cosh(Real(z)) */ feclearexcept(FE_UNDERFLOW); /* underflows don't matter */ feupdateenv(&env); /* update environment */ return w; } /**************************************************************************** double complex ccosh(double complex z) returns the complex hyperbolic cosine of its argument. The algorithm is based upon the identity: cosh(x + i*y) = cos(y)*cosh(x) + i*sin(y)*sinh(x). Signaling of spurious overflows, underflows, and invalids is avoided where possible. Calls: sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex ccosh ( double complex z ) { fenv_t env; int ifoo; double sinval, cosval; double complex w; ifoo = feholdexcept(&env); /* save environment, clear flags */ sinval = sin(Imag(z)); /* sine of imaginary part of argument */ cosval = cos(Imag(z)); /* cosine of imaginary part of argument */ Real(w) = coshmul(Real(z),cosval); /* real result = cosval*cosh(Real(z)) */ Imag(w) = sinhmul(Real(z),sinval); /* imag result = sinval*sinh(Real(z)) */ feclearexcept(FE_UNDERFLOW); /* underflows don't matter */ feupdateenv(&env); /* update environment */ return w; } /**************************************************************************** double complex cexp(double complex z) returns the complex exponential of its argument. The algorithm is based upon the identity: exp(x + i*y) = cos(y)*exp(x) + i*sin(y)*exp(x). Signaling of spurious overflows and invalids is avoided where possible. CONSTANTS: FPKEXP709 = exp(709.0) to double precision Calls: sin, cos, exp, fpclassify, and scalb. Called by: cepwry, cxpwri, cxpwre, and cxpwry ****************************************************************************/ double complex cexp ( double complex z ) { double sinval, cosval, expval, exparg; double complex w; sinval = sin(Imag(z)); /* sine of imaginary part of argument */ cosval = cos(Imag(z)); /* cosine of imaginary part of argument */ if (Real(z) <= 709.0) { /* exp(Real(z)) is finite */ expval = exp(Real(z)); /* evaluate exponential */ Real(w) = cosval*expval; /* real result = cos(Imag(z))*exp(Real(z)) */ Imag(w) = sinval*expval; /* imag result = sin(Imag(z))*exp(Real(z)) */ } else if (fpclassify(Real(z)) < FP_ZERO) { /* Real(z) = +INF or a NaN */ Real(w) = cosval*Real(z); /* deserved invalid may occur */ Imag(w) = sinval*Real(z); /* deserved invalid may occur */ } else if (Real(z) > 1460.0) { /* probable overflow case */ Real(w) = scalb(cosval,2100); Imag(w) = scalb(sinval,2100); } else { /* exp(Real(z)) overflows but product with sin or cos may not */ Real(w) = cosval*FPKEXP709.dval; /* initialize real result */ Imag(w) = sinval*FPKEXP709.dval; /* initialize imag result */ exparg = Real(z) - 709.0; /* initialize reduced exponent argument */ while (exparg > 709.0) { /* exponential reduction loop */ Real(w) *= FPKEXP709.dval; Imag(w) *= FPKEXP709.dval; exparg -= 709.0; } expval = exp(exparg); /* final exponential value */ Real(w) *= expval; /* final multiplication steps */ Imag(w) *= expval; } return w; /* done */ } /**************************************************************************** double complex cpow(double complex x, double complex y) returns the complex result of x^y. The algorithm is based upon the identity: x^y = cexp(y*clog(x)). Calls: clog, cexp. ****************************************************************************/ double complex cpow ( double complex x, double complex y ) /* (complex x)^(complex y) */ { double complex logval,z; logval = clog(x); /* complex logarithm of x */ Real(z) = Real(y)*Real(logval) - Imag(y)*Imag(logval); /* multiply by y */ Imag(z) = Real(y)*Imag(logval) + Imag(y)*Real(logval); return (cexp(z)); /* return complex exponential */ } /**************************************************************************** double complex ctanh(double complex x) returns the complex hyperbolic tangent of its argument. The algorithm is from Kahan's paper and is based on the identity: tanh(x+i*y) = (sinh(2*x) + i*sin(2*y))/(cosh(2*x) + cos(2*y)) = (cosh(x)*sinh(x)*cscsq + i*tan(y))/(1+cscsq*sinh(x)*sinh(x)), where cscsq = 1/(cos(y)*cos(y). For large values of ze.re, spurious overflow and invalid signaling is avoided. CONSTANTS: FPKASINHOM4 = asinh(nextafterd(+infinity,0.0))/4.0 to double precision FPKINF = +infinity Calls: tan, sinh, sqrt, fabs, feholdexcept, feraiseexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex ctanh( double complex z ) { fenv_t env; int ifoo; double tanval, beta, sinhval, coshval, denom; double complex w; ifoo = feholdexcept(&env); /* save environment and clear flags */ if (fabs(Real(z)) > FPKASINHOM4.dval) { /* avoid overflow for large |Real(z)| */ Real(w) = copysign(1.0,Real(z)); /* real result has unit magnitude */ Imag(w) = copysign(0.0,Imag(z)); /* imag result is signed zero */ if (fabs(Real(z)) != FPKINF.dval) /* set inexact for finite Real(z) */ feraiseexcept(FE_INEXACT); feupdateenv(&env); /* update environment */ } /* end large |Real(z)| case */ else { /* usual case */ tanval = tan(Imag(z)); /* evaluate tangent */ feclearexcept(FE_DIVBYZERO); /* in case tangent is infinite */ feupdateenv(&env); /* update environment */ beta = 1.0 + tanval*tanval; /* 1/(cos(Imag(z)))^2 */ sinhval = sinh(Real(z)); /* evaluate sinh */ coshval = sqrt(1.0+sinhval*sinhval); /* evaluate cosh */ if (fabs(tanval) == FPKINF.dval) { /* infinite tangent */ Real(w) = coshval/sinhval; Imag(w) = 1.0/tanval; } else { /* finite tangent */ denom = 1.0 + beta*sinhval*sinhval; Real(w) = beta*coshval*sinhval/denom; Imag(w) = tanval/denom; } } /* end usual case */ return w; /* done */ } /**************************************************************************** double complex ctan(double complex x) returns the complex hyperbolic tangent of its argument. The algorithm is from Kahan's paper and is based on the identity: tan(x + i*y) = (sin(2*x) + i*sinh(2*y))/(cos(2*x) + cosh(2*y)) = (tan(x)+i*cosh(y)*sinh(y)*cscsq)/(1+cscsq*sinh(y)*sinh(y)), where cscsq = 1/(cos(x)*cos(x). For large values of ze.im, spurious overflow and invalid signaling is avoided. CONSTANTS: FPKASINHOM4 = asinh(nextafterd(+infinity,0.0))/4.0 to double precision FPKINF = +infinity Calls: tan, sinh, sqrt, fabs, feholdexcept, feraiseexcept, feclearexcept, and feupdateenv. ****************************************************************************/ double complex ctan( double complex z ) { fenv_t env; int ifoo; double tanval, beta, sinhval, coshval, denom; double complex w; ifoo = feholdexcept(&env); /* save environment and clear flags */ if (fabs(Imag(z)) > FPKASINHOM4.dval) { /* avoid overflow for large |Imag(z)| */ Real(w) = copysign(0.0,Real(z)); /* real result has zero magnitude */ Imag(w) = copysign(1.0,Imag(z)); /* imag result has unit magnitude */ if (fabs(Imag(z)) != FPKINF.dval) /* set inexact for finite Real(z) */ feraiseexcept(FE_INEXACT); feupdateenv(&env); /* update environment */ } /* end large |Real(z)| case */ else { /* usual case */ tanval = tan(Real(z)); /* evaluate tangent */ feclearexcept(FE_DIVBYZERO); /* in case tangent is infinite */ feupdateenv(&env); /* update environment */ beta = 1.0 + tanval*tanval; /* 1/(cos(Real(z)))^2 */ sinhval = sinh(Imag(z)); /* evaluate sinh */ coshval = sqrt(1.0+sinhval*sinhval); /* evaluate cosh */ if (fabs(tanval) == FPKINF.dval) { /* infinite tangent */ Real(w) = 1.0/tanval; Imag(w) = coshval/sinhval; } else { /* finite tangent */ denom = 1.0 + beta*sinhval*sinhval; Real(w) = tanval/denom; Imag(w) = beta*coshval*sinhval/denom; } } /* end usual case */ return w; /* done */ } /**************************************************************************** static double arcsinh(double x) returns the real inverse hyperbolic sine of its real argument. The implementation is based on the identity: arcsinh(x) = log(x + sqrt(x*x + 1)) = log1p(x + x/(1/x + sqrt(1 + 1/(x*x)))) CONSTANTS: FPK2M33 = 2.0^(-33) FPK2P54 = 2.0^54 FPKLOG2 = log(2.0) to double precision FPKTINY = 2.0^(-1022) Calls: fpclassify, fabs, log, log1p, sqrt, copysign, feraiseexcept. Called by: casinh, casin, cacosh, and cacos. ****************************************************************************/ static double arcsinh ( double x ) { double y; int iclass; iclass = fpclassify(x); /* classify argument */ if (iclass == FP_NORMAL) { y = fabs(x); /* normal argument */ if (y < FPK2M33.dval) /* tiny normal argument */ feraiseexcept(FE_INEXACT); /* raise inexact and return arg */ else if (y > FPK2P54.dval) /* huge magnitude argument */ y = log(y) + FPKLOG2.dval; /* prevents spurious overflows */ else /* usual case */ y = log1p(y + y/(1.0/y + sqrt(1.0 + 1.0/(y*y)))); return copysign(y,x); /* fix sign */ } /* end normal case */ else if (iclass < FP_NORMAL) /* zero, INF, or NaN argument */ return (x + x); /* propagate argument */ else { /* subnormal argument */ feraiseexcept(FE_UNDERFLOW + FE_INEXACT); /* signal underflow/inex */ return (x); /* return argument */ } } /**************************************************************************** double complex casin(double complex z) returns the complex inverse sine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(casin(z)) = atan (real(z)/real(csqrt(1.0-z)*csqrt(1.0+z))) imag(casin(z)) = arcsinh(imag(csqrt(1.0-cconj(z))*csqrt(1.0+z))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv. ****************************************************************************/ double complex casin ( double complex z ) { double complex zp1, zm, zm1, w; fenv_t env; int ifoo; Real(zp1) = 1.0 + Real(z); /* evaluate zp1 = csqrt(1.0+z) */ Imag(zp1) = Imag(z); zp1 = csqrt(zp1); Real(zm) = 1.0 - Real(z); /* evaluate zm1 = csqrt(1.0-z) */ Imag(zm) = -Imag(z); zm1 = csqrt(zm); ifoo = feholdexcept(&env); /* save environ., clr flags/enables */ Real(w) = atan(Real(z)/(Real(zp1)*Real(zm1) -Imag(zp1)*Imag(zm1))); /* real result */ feclearexcept(FE_DIVBYZERO); /* clear spurious DIVBYZERO exception */ feupdateenv(&env); /* restore environment with new flags */ Imag(zm) = Imag(z); /* evaluate zm1 = csqrt(1.0-cconj(z)) */ zm1 = csqrt(zm); Imag(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1)); /* imag result */ return w; } /**************************************************************************** double complex casinh(double complex z) returns the complex inverse hyperbolic sine of its argument. The algorithm is from Kahan's paper and is based on the formula: casinh(z) = -i*casin(i*z). Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv. ****************************************************************************/ double complex casinh ( double complex z ) { double complex zp1, zm, zm1, w; fenv_t env; int ifoo; Real(zp1) = 1.0 - Imag(z); /* evaluate zp1 = csqrt(1.0+i*z) */ Imag(zp1) = Real(z); zp1 = csqrt(zp1); Real(zm) = 1.0 + Imag(z); /* evaluate zm1 = csqrt(1.0-cconj(i*z)) */ Imag(zm) = Real(z); zm1 = csqrt(zm); Real(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1)); /* imag result */ Imag(zm) = -Real(z); /* evaluate zm1 = csqrt(1.0-cconj(z)) */ zm1 = csqrt(zm); ifoo = feholdexcept(&env); /* save environ., clr flags/enables */ Imag(w) = atan(Imag(z)/(Real(zp1)*Real(zm1) -Imag(zp1)*Imag(zm1))); /* real result */ feclearexcept(FE_DIVBYZERO); /* clear spurious DIVBYZERO exception */ feupdateenv(&env); /* restore environment with new flags */ return w; } /**************************************************************************** double complex cacos(double complex z) returns the complex inverse cosine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(cacos(z)) = 2.0*atan(real(csqrt(1.0-z)/real(csqrt(1.0+z)))) imag(cacos(z)) = arcsinh(imag(csqrt(1.0-z)*csqrt(cconj(1.0+z)))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv. ****************************************************************************/ double complex cacos ( double complex z ) { double complex zp, zp1, zm1, w; fenv_t env; int ifoo; Real(zp) = 1.0 + Real(z); /* zp1 = csqrt(1.0 + z) */ Imag(zp) = Imag(z); zp1 = csqrt(zp); Real(zm1) = 1.0 - Real(z); /* zm1 = csqrt(1.0 - z) */ Imag(zm1) = -Imag(z); zm1 = csqrt(zm1); ifoo = feholdexcept(&env); /* save environment, clr flags/enables */ Real(w) = 2.0*atan(Real(zm1)/Real(zp1)); /* real result */ feclearexcept(FE_DIVBYZERO); /* clr possible spurious div-by-zero flag */ feupdateenv(&env); /* update environment with new flags */ Imag(zp) = -Imag(z); /* zp1 = csqrt(1.0 + cconj(z)) */ zp1 = csqrt(zp); Imag(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1)); /* imag result */ return w; } /**************************************************************************** double complex cacosh(double complex z) returns the complex inverse hyperbolic`cosine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(cacosh(z)) = arcsinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0))) imag(cacosh(z)) = 2.0*atan(imag(csqrt(z-1.0)/real(csqrt(z+1.0)))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv. ****************************************************************************/ double complex cacosh ( double complex z ) { double complex zp1, zm, zm1, w; fenv_t env; int ifoo; Real(zp1) = Real(z) + 1.0; /* zp1 = csqrt(z + 1.0) */ Imag(zp1) = Imag(z); zp1 = csqrt(zp1); Real(zm) = Real(z) - 1.0; /* zm1 = csqrt(z - 1.0) */ Imag(zm) = Imag(z); zm1 = csqrt(zm); ifoo = feholdexcept(&env); /* save environment, clr flags/enables */ Imag(w) = 2.0*atan(Imag(zm1)/Real(zp1)); /* imag result */ feclearexcept(FE_DIVBYZERO); /* clr possible spurious div-by-zero flag */ feupdateenv(&env); /* restore environment with new flags */ Imag(zm) = -Imag(z); /* zm1 = csqrt(cconj(z) - 1.0) */ zm1 = csqrt(zm); Real(w) = arcsinh(Real(zp1)*Real(zm1) - Imag(zp1)*Imag(zm1)); /* real result */ return w; } /**************************************************************************** double complex catan(double complex z) returns the complex inverse tangent of its argument. The algorithm is from Kahan's paper and is based on the formula: catan(z) = i*(clog(1.0-i*z) - clog(1+i*z))/2.0. CONSTANTS: FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0 FPKRHO = 1.0/FPKTHETA FPKPI2 = pi/2.0 FPKPI = pi Calls: copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg. ****************************************************************************/ double complex catan ( double complex z ) { double complex ctemp, w; double t1, t2, xi, eta, beta; xi = -Imag(z); beta = copysign(1.0,xi); /* copes with unsigned zero */ Imag(z) = -beta*Real(z); /* transform real & imag components */ Real(z) = beta*xi; if ((Real(z) > FPKTHETA.dval) || (fabs(Imag(z)) > FPKTHETA.dval)) { xi = copysign(FPKPI2.dval,Imag(z)); /* avoid spurious overflow */ ctemp = xdivc(1.0,z); eta = Real(ctemp); } else if (Real(z) == 1.0) { t1 = fabs(Imag(z)) + FPKRHO.dval; xi = log(sqrt(sqrt(4.0 + t1*t1))/sqrt(fabs(Imag(z)))); eta = 0.5*copysign(FPKPI.dval-atan(2.0/(fabs(Imag(z))+FPKRHO.dval)),Imag(z)); } else { /* usual case */ t2 = fabs(Imag(z)) + FPKRHO.dval; t1 = 1.0 - Real(z); t2 = t2*t2; xi = 0.25*log1p(4.0*Real(z)/(t1*t1 + t2)); Real(ctemp) = (1.0 - Real(z))*(1.0 + Real(z)) - t2; Imag(ctemp) = Imag(z) + Imag(z); eta = 0.5*carg(ctemp); } Real(w) = -beta*eta; /* fix up signs of result */ Imag(w) = -beta*xi; return w; } /**************************************************************************** double complex catanh(double complex z) returns the complex inverse hyperbolic tangent of its argument. The algorithm is from Kahan's paper and is based on the formula: catan(z) = (clog(1.0 + z) - clog(1 - z))/2.0. CONSTANTS: FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0 FPKRHO = 1.0/FPKTHETA FPKPI2 = pi/2.0 FPKPI = pi Calls: copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg. ****************************************************************************/ double complex catanh( double complex z ) { double complex ctemp, w; double t1, t2, xi, eta, beta; beta = copysign(1.0,Real(z)); /* copes with unsigned zero */ Imag(z) = -beta*Imag(z); /* transform real & imag components */ Real(z) = beta*Real(z); if ((Real(z) > FPKTHETA.dval) || (fabs(Imag(z)) > FPKTHETA.dval)) { eta = copysign(FPKPI2.dval,Imag(z)); /* avoid overflow */ ctemp = xdivc(1.0,z); xi = Real(ctemp); } else if (Real(z) == 1.0) { t1 = fabs(Imag(z)) + FPKRHO.dval; xi = log(sqrt(sqrt(4.0 + t1*t1))/sqrt(fabs(Imag(z)))); eta = 0.5*copysign(FPKPI.dval-atan(2.0/(fabs(Imag(z))+FPKRHO.dval)),Imag(z)); } else { /* usual case */ t2 = fabs(Imag(z)) + FPKRHO.dval; t1 = 1.0 - Real(z); t2 = t2*t2; xi = 0.25*log1p(4.0*Real(z)/(t1*t1 + t2)); Real(ctemp) = (1.0 - Real(z))*(1.0 + Real(z)) - t2; Imag(ctemp) = Imag(z) + Imag(z); eta = 0.5*carg(ctemp); } Real(w) = beta*xi; /* fix up signs of result */ Imag(w) = -beta*eta; return w; } /* conj(), creal(), and cimag() are gcc built ins. */ double complex cproj( double complex z ) { double complex S; double zR = Real(z), zI = Imag(z), scale; if (isnan(zR) && isnan(zI)) return z; if (isinf(zR) || isinf(zI)) { Real(S) = INFINITY; Imag(S) = copysign(0.0, zI); return S; } scale = 1.0 + zR * zR + zI * zI; Real(S) = (zR + zR)/scale; Imag(S) = (zI + zI)/scale; return S; } float complex cacosf( float complex z ){ return (float complex) cacos((double complex) z); } float complex casinf( float complex z ){ return (float complex) casin((double complex) z); } float complex catanf( float complex z ){ return (float complex) catan((double complex) z); } float complex ccosf( float complex z ){ return (float complex) ccos((double complex) z); } float complex csinf( float complex z ){ return (float complex) csin((double complex) z); } float complex ctanf( float complex z ){ return (float complex) ctan((double complex) z); } float complex cacoshf( float complex z ){ return (float complex) cacosh((double complex) z); } float complex casinhf( float complex z ){ return (float complex) casinh((double complex) z); } float complex catanhf( float complex z ){ return (float complex) catanh((double complex) z); } float complex ccoshf( float complex z ){ return (float complex) ccosh((double complex) z); } float complex csinhf( float complex z ){ return (float complex) csinh((double complex) z); } float complex ctanhf( float complex z ){ return (float complex) ctanh((double complex) z); } float complex cexpf( float complex z ){ return (float complex) cexp((double complex) z); } float complex clogf( float complex z ){ return (float complex) clog((double complex) z); } float complex cpowf( float complex x, float complex y ){ return (float complex) cpow((double complex) x, (double complex) y); } float complex csqrtf( float complex z ){ return (float complex) csqrt((double complex) z); } float cargf( float complex z ){ return (float) carg((double complex) z); } float complex cprojf( float complex z ){ return (float complex) cproj((double complex) z); }