Стандарт языка Си С99 TC (1113411), страница 80
Текст из файла (страница 80)
Notethat the imaginary unit I has imaginary type (see G.6).#include <math.h>#include <complex.h>/* Multiply z * w ... */double complex _Cmultd(double complex z, double complex w){#pragma STDC FP_CONTRACT OFFdouble a, b, c, d, ac, bd, ad, bc, x, y;a = creal(z); b = cimag(z);c = creal(w); d = cimag(w);ac = a * c;bd = b * d;ad = a * d;bc = b * c;x = ac - bd; y = ad + bc;if (isnan(x) && isnan(y)) {/* Recover infinities that computed as NaN+iNaN ...
*/int recalc = 0;if ( isinf(a) || isinf(b) ) { // z is infinite/* "Box" the infinity and change NaNs in the other factor to 0 */a = copysign(isinf(a) ? 1.0 : 0.0, a);b = copysign(isinf(b) ? 1.0 : 0.0, b);if (isnan(c)) c = copysign(0.0, c);if (isnan(d)) d = copysign(0.0, d);recalc = 1;}if ( isinf(c) || isinf(d) ) { // w is infinite/* "Box" the infinity and change NaNs in the other factor to 0 */c = copysign(isinf(c) ? 1.0 : 0.0, c);d = copysign(isinf(d) ? 1.0 : 0.0, d);if (isnan(a)) a = copysign(0.0, a);if (isnan(b)) b = copysign(0.0, b);recalc = 1;}if (!recalc && (isinf(ac) || isinf(bd) ||isinf(ad) || isinf(bc))) {/* Recover infinities from overflow by changing NaNs to 0 ...
*/if (isnan(a)) a = copysign(0.0, a);if (isnan(b)) b = copysign(0.0, b);if (isnan(c)) c = copysign(0.0, c);if (isnan(d)) d = copysign(0.0, d);recalc = 1;}if (recalc) {470IEC 60559-compatible complex arithmetic§G.5.1WG14/N1256Committee Draft — Septermber 7, 2007ISO/IEC 9899:TC3x = INFINITY * ( a * c - b * d );y = INFINITY * ( a * d + b * c );}}return x + I * y;}7This implementation achieves the required treatment of infinities at the cost of only one isnan test inordinary (finite) cases. It is less than ideal in that undue overflow and underflow may occur.8EXAMPLE 2Division of two double _Complex operands could be implemented as follows.#include <math.h>#include <complex.h>/* Divide z / w ... */double complex _Cdivd(double complex z, double complex w){#pragma STDC FP_CONTRACT OFFdouble a, b, c, d, logbw, denom, x, y;int ilogbw = 0;a = creal(z); b = cimag(z);c = creal(w); d = cimag(w);logbw = logb(fmax(fabs(c), fabs(d)));if (isfinite(logbw)) {ilogbw = (int)logbw;c = scalbn(c, -ilogbw); d = scalbn(d, -ilogbw);}denom = c * c + d * d;x = scalbn((a * c + b * d) / denom, -ilogbw);y = scalbn((b * c - a * d) / denom, -ilogbw);/* Recover infinities and zeros that computed as NaN+iNaN;*//* the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
*/if (isnan(x) && isnan(y)) {if ((denom == 0.0) &&(!isnan(a) || !isnan(b))) {x = copysign(INFINITY, c) * a;y = copysign(INFINITY, c) * b;}else if ((isinf(a) || isinf(b)) &&isfinite(c) && isfinite(d)) {a = copysign(isinf(a) ? 1.0 : 0.0,b = copysign(isinf(b) ? 1.0 : 0.0,x = INFINITY * ( a * c + b * d );y = INFINITY * ( b * c - a * d );}else if (isinf(logbw) &&isfinite(a) && isfinite(b)) {c = copysign(isinf(c) ? 1.0 : 0.0,d = copysign(isinf(d) ? 1.0 : 0.0,x = 0.0 * ( a * c + b * d );y = 0.0 * ( b * c - a * d );§G.5.1IEC 60559-compatible complex arithmetica);b);c);d);471ISO/IEC 9899:TC3Committee Draft — Septermber 7, 2007WG14/N1256}}return x + I * y;}9Scaling the denominator alleviates the main overflow and underflow problem, which is more serious thanfor multiplication.
In the spirit of the multiplication example above, this code does not defend againstoverflow and underflow in the calculation of the numerator. Scaling with the scalbn function, instead ofwith division, provides better roundoff characteristics.G.5.2 Additive operatorsSemantics1If both operands have imaginary type, then the result has imaginary type. (If one operandhas real type and the other operand has imaginary type, or if either operand has complextype, then the result has complex type.)2In all cases the result and floating-point exception behavior of a + or - operator is definedby the usual mathematical formula:uivu + ivxx±ux ± iv(x ± u) ± iviy±u + iyi(y ± v)±u + i(y ± v)(x ± u) + iyx + i(y ± v)(x ± u) + i(y ± v)+ or −x + iyG.6 Complex arithmetic <complex.h>1The macrosimaginaryand_Imaginary_Iare defined, respectively, as _Imaginary and a constant expression of type constfloat _Imaginary with the value of the imaginary unit.
The macroIis defined to be _Imaginary_I (not _Complex_I as stated in 7.3). Notwithstandingthe provisions of 7.1.3, a program may undefine and then perhaps redefine the macroimaginary.2This subclause contains specifications for the <complex.h> functions that areparticularly suited to IEC 60559 implementations.
For families of functions, thespecifications apply to all of the functions even though only the principal function is472IEC 60559-compatible complex arithmetic§G.6WG14/N1256Committee Draft — Septermber 7, 2007ISO/IEC 9899:TC3shown. Unless otherwise specified, where the symbol ‘‘±’’ occurs in both an argumentand the result, the result has the same sign as the argument.3The functions are continuous onto both sides of their branch cuts, taking into account thesign of zero. For example, csqrt(−2 ± i0) = ±i√ 2.4Since complex and imaginary values are composed of real values, each function may beregarded as computing real values from real values.
Except as noted, the functions treatreal infinities, NaNs, signed zeros, subnormals, and the floating-point exception flags in amanner consistent with the specifications for real functions in F.9.326)5The functions cimag, conj, cproj, and creal are fully specified for allimplementations, including IEC 60559 ones, in 7.3.9.
These functions raise no floatingpoint exceptions.6Each of the functions cabs and carg is specified by a formula in terms of a realfunction (whose special cases are covered in annex F):cabs(x + iy) = hypot(x, y)carg(x + iy) = atan2(y, x)7Each of the functions casin, catan, ccos, csin, and ctan is specified implicitly bya formula in terms of other complex functions (whose special cases are specified below):casin(z)catan(z)ccos(z)csin(z)ctan(z)=====−i casinh(iz)−i catanh(iz)ccosh(iz)−i csinh(iz)−i ctanh(iz)8For the other functions, the following subclauses specify behavior for special cases,including treatment of the ‘‘invalid’’ and ‘‘divide-by-zero’’ floating-point exceptions.
Forfamilies of functions, the specifications apply to all of the functions even though only theprincipal function is shown. For a function f satisfying f (conj(z)) = conj( f (z)), thespecifications for the upper half-plane imply the specifications for the lower half-plane; ifthe function f is also either even, f (−z) = f (z), or odd, f (−z) = − f (z), then thespecifications for the first quadrant imply the specifications for the other three quadrants.9In the following subclauses, cis(y) is defined as cos(y) + i sin(y).326) As noted in G.3, a complex value with at least one infinite part is regarded as an infinity even if itsother part is a NaN.§G.6IEC 60559-compatible complex arithmetic473ISO/IEC 9899:TC3Committee Draft — Septermber 7, 2007WG14/N1256G.6.1 Trigonometric functionsG.6.1.1 The cacos functions1— cacos(conj(z)) = conj(cacos(z)).— cacos(±0 + i0) returns π /2 − i0.— cacos(±0 + iNaN) returns π /2 + iNaN.— cacos(x + i ∞) returns π /2 − i ∞, for finite x.— cacos(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for nonzero finite x.— cacos(−∞ + iy) returns π − i ∞, for positive-signed finite y.— cacos(+∞ + iy) returns +0 − i ∞, for positive-signed finite y.— cacos(−∞ + i ∞) returns 3π /4 − i ∞.— cacos(+∞ + i ∞) returns π /4 − i ∞.— cacos(±∞ + iNaN) returns NaN ± i ∞ (where the sign of the imaginary part of theresult is unspecified).— cacos(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite y.— cacos(NaN + i ∞) returns NaN − i ∞.— cacos(NaN + iNaN) returns NaN + iNaN.G.6.2 Hyperbolic functionsG.6.2.1 The cacosh functions1— cacosh(conj(z)) = conj(cacosh(z)).— cacosh(±0 + i0) returns +0 + iπ /2.— cacosh(x + i ∞) returns +∞ + iπ /2, for finite x.— cacosh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for finite x.— cacosh(−∞ + iy) returns +∞ + iπ , for positive-signed finite y.— cacosh(+∞ + iy) returns +∞ + i0, for positive-signed finite y.— cacosh(−∞ + i ∞) returns +∞ + i3π /4.— cacosh(+∞ + i ∞) returns +∞ + iπ /4.— cacosh(±∞ + iNaN) returns +∞ + iNaN.474IEC 60559-compatible complex arithmetic§G.6.2.1WG14/N1256Committee Draft — Septermber 7, 2007ISO/IEC 9899:TC3— cacosh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for finite y.— cacosh(NaN + i ∞) returns +∞ + iNaN.— cacosh(NaN + iNaN) returns NaN + iNaN.G.6.2.2 The casinh functions1— casinh(conj(z)) = conj(casinh(z)) and casinh is odd.— casinh(+0 + i0) returns 0 + i0.— casinh(x + i ∞) returns +∞ + iπ /2 for positive-signed finite x.— casinh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for finite x.— casinh(+∞ + iy) returns +∞ + i0 for positive-signed finite y.— casinh(+∞ + i ∞) returns +∞ + iπ /4.— casinh(+∞ + iNaN) returns +∞ + iNaN.— casinh(NaN + i0) returns NaN + i0.— casinh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for finite nonzero y.— casinh(NaN + i ∞) returns ±∞ + iNaN (where the sign of the real part of the resultis unspecified).— casinh(NaN + iNaN) returns NaN + iNaN.G.6.2.3 The catanh functions1— catanh(conj(z)) = conj(catanh(z)) and catanh is odd.— catanh(+0 + i0) returns +0 + i0.— catanh(+0 + iNaN) returns +0 + iNaN.— catanh(+1 + i0) returns +∞ + i0 and raises the ‘‘divide-by-zero’’ floating-pointexception.— catanh(x + i ∞) returns +0 + iπ /2, for finite positive-signed x.— catanh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for nonzero finite x.— catanh(+∞ + iy) returns +0 + iπ /2, for finite positive-signed y.— catanh(+∞ + i ∞) returns +0 + iπ /2.— catanh(+∞ + iNaN) returns +0 + iNaN.§G.6.2.3IEC 60559-compatible complex arithmetic475ISO/IEC 9899:TC3Committee Draft — Septermber 7, 2007WG14/N1256— catanh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’floating-point exception, for finite y.— catanh(NaN + i ∞) returns ±0 + iπ /2 (where the sign of the real part of the result isunspecified).— catanh(NaN + iNaN) returns NaN + iNaN.G.6.2.4 The ccosh functions1— ccosh(conj(z)) = conj(ccosh(z)) and ccosh is even.— ccosh(+0 + i0) returns 1 + i0.— ccosh(+0 + i ∞) returns NaN ± i0 (where the sign of the imaginary part of theresult is unspecified) and raises the ‘‘invalid’’ floating-point exception.— ccosh(+0 + iNaN) returns NaN ± i0 (where the sign of the imaginary part of theresult is unspecified).— ccosh(x + i ∞) returns NaN + iNaN and raises the ‘‘invalid’’ floating-pointexception, for finite nonzero x.— ccosh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite nonzero x.— ccosh(+∞ + i0) returns +∞ + i0.— ccosh(+∞ + iy) returns +∞ cis(y), for finite nonzero y.— ccosh(+∞ + i ∞) returns ±∞ + iNaN (where the sign of the real part of the result isunspecified) and raises the ‘‘invalid’’ floating-point exception.— ccosh(+∞ + iNaN) returns +∞ + iNaN.— ccosh(NaN + i0) returns NaN ± i0 (where the sign of the imaginary part of theresult is unspecified).— ccosh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for all nonzero numbers y.— ccosh(NaN + iNaN) returns NaN + iNaN.G.6.2.5 The csinh functions1— csinh(conj(z)) = conj(csinh(z)) and csinh is odd.— csinh(+0 + i0) returns +0 + i0.— csinh(+0 + i ∞) returns ±0 + iNaN (where the sign of the real part of the result isunspecified) and raises the ‘‘invalid’’ floating-point exception.— csinh(+0 + iNaN) returns ±0 + iNaN (where the sign of the real part of the result isunspecified).476IEC 60559-compatible complex arithmetic§G.6.2.5WG14/N1256Committee Draft — Septermber 7, 2007ISO/IEC 9899:TC3— csinh(x + i ∞) returns NaN + iNaN and raises the ‘‘invalid’’ floating-pointexception, for positive finite x.— csinh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite nonzero x.— csinh(+∞ + i0) returns +∞ + i0.— csinh(+∞ + iy) returns +∞ cis(y), for positive finite y.— csinh(+∞ + i ∞) returns ±∞ + iNaN (where the sign of the real part of the result isunspecified) and raises the ‘‘invalid’’ floating-point exception.— csinh(+∞ + iNaN) returns ±∞ + iNaN (where the sign of the real part of the resultis unspecified).— csinh(NaN + i0) returns NaN + i0.— csinh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for all nonzero numbers y.— csinh(NaN + iNaN) returns NaN + iNaN.G.6.2.6 The ctanh functions1— ctanh(conj(z)) = conj(ctanh(z))and ctanh is odd.— ctanh(+0 + i0) returns +0 + i0.— ctanh(x + i ∞) returns NaN + iNaN and raises the ‘‘invalid’’ floating-pointexception, for finite x.— ctanh(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite x.— ctanh(+∞ + iy) returns 1 + i0 sin(2y), for positive-signed finite y.— ctanh(+∞ + i ∞) returns 1 ± i0 (where the sign of the imaginary part of the resultis unspecified).— ctanh(+∞ + iNaN) returns 1 ± i0 (where the sign of the imaginary part of theresult is unspecified).— ctanh(NaN + i0) returns NaN + i0.— ctanh(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for all nonzero numbers y.— ctanh(NaN + iNaN) returns NaN + iNaN.§G.6.2.6IEC 60559-compatible complex arithmetic477ISO/IEC 9899:TC3Committee Draft — Septermber 7, 2007WG14/N1256G.6.3 Exponential and logarithmic functionsG.6.3.1 The cexp functions1— cexp(conj(z)) = conj(cexp(z)).— cexp(±0 + i0) returns 1 + i0.— cexp(x + i ∞) returns NaN + iNaN and raises the ‘‘invalid’’ floating-pointexception, for finite x.— cexp(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite x.— cexp(+∞ + i0) returns +∞ + i0.— cexp(−∞ + iy) returns +0 cis(y), for finite y.— cexp(+∞ + iy) returns +∞ cis(y), for finite nonzero y.— cexp(−∞ + i ∞) returns ±0 ± i0 (where the signs of the real and imaginary parts ofthe result are unspecified).— cexp(+∞ + i ∞) returns ±∞ + iNaN and raises the ‘‘invalid’’ floating-pointexception (where the sign of the real part of the result is unspecified).— cexp(−∞ + iNaN) returns ±0 ± i0 (where the signs of the real and imaginary partsof the result are unspecified).— cexp(+∞ + iNaN) returns ±∞ + iNaN (where the sign of the real part of the resultis unspecified).— cexp(NaN + i0) returns NaN + i0.— cexp(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for all nonzero numbers y.— cexp(NaN + iNaN) returns NaN + iNaN.G.6.3.2 The clog functions1— clog(conj(z)) = conj(clog(z)).— clog(−0 + i0) returns −∞ + iπ and raises the ‘‘divide-by-zero’’ floating-pointexception.— clog(+0 + i0) returns −∞ + i0 and raises the ‘‘divide-by-zero’’ floating-pointexception.— clog(x + i ∞) returns +∞ + iπ /2, for finite x.— clog(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite x.478IEC 60559-compatible complex arithmetic§G.6.3.2WG14/N1256Committee Draft — Septermber 7, 2007ISO/IEC 9899:TC3— clog(−∞ + iy) returns +∞ + iπ , for finite positive-signed y.— clog(+∞ + iy) returns +∞ + i0, for finite positive-signed y.— clog(−∞ + i ∞) returns +∞ + i3π /4.— clog(+∞ + i ∞) returns +∞ + iπ /4.— clog(±∞ + iNaN) returns +∞ + iNaN.— clog(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite y.— clog(NaN + i ∞) returns +∞ + iNaN.— clog(NaN + iNaN) returns NaN + iNaN.G.6.4 Power and absolute-value functionsG.6.4.1 The cpow functions1The cpow functions raise floating-point exceptions if appropriate for the calculation ofthe parts of the result, and may raise spurious exceptions.327)G.6.4.2 The csqrt functions1— csqrt(conj(z)) = conj(csqrt(z)).— csqrt(±0 + i0) returns +0 + i0.— csqrt(x + i ∞) returns +∞ + i ∞, for all x (including NaN).— csqrt(x + iNaN) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite x.— csqrt(−∞ + iy) returns +0 + i ∞, for finite positive-signed y.— csqrt(+∞ + iy) returns +∞ + i0, for finite positive-signed y.— csqrt(−∞ + iNaN) returns NaN ± i ∞ (where the sign of the imaginary part of theresult is unspecified).— csqrt(+∞ + iNaN) returns +∞ + iNaN.— csqrt(NaN + iy) returns NaN + iNaN and optionally raises the ‘‘invalid’’ floatingpoint exception, for finite y.— csqrt(NaN + iNaN) returns NaN + iNaN.327) This allows cpow( z , c ) to be implemented as cexp(cimplementations that treat special cases more carefully.§G.6.4.2clog( z )) without precludingIEC 60559-compatible complex arithmetic479ISO/IEC 9899:TC3Committee Draft — Septermber 7, 2007WG14/N1256G.7 Type-generic math <tgmath.h>1Type-generic macros that accept complex arguments also accept imaginary arguments.