/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:23 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "zwofz.h" #include #include /* PARAMETER translations */ #define LN2 0.6931471805599453094172321214581765680755e0 #define RSQPD2 1.128379167095512573896158903121545171688e0 /* end of PARAMETER translations */ void /*FUNCTION*/ zwofz( double z[], double w[], long *flag) { LOGICAL32 a; long int _l0, i, j, kapn, n, nu; double c, h, h2, qlamda, qrho, rx, ry, sx, sy, tx, ty, u, u1, u2, uv, v, v1, v2, w1, x, xabs, xquad, xsum, y, yabs, yquad, ysum; static double inogxm, lgund, mxgoni, sqlgov; static double lgov2 = -1.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const W = &w[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>>2002-11-06 ZWOFZ Krogh Corrected comments. *>>1996-03-30 ZWOFZ Krogh Added external stmt., removed INT in type st. *>>1992-10-13 ZWOFZ WVS Improve efficiency and avoid underflow. *>>1992-03-13 ZWOFZ FTK Removed implicit statements. *>>1991-08-23 ZWOFZ WV Snyder Initial adaptation to Math77 * * Algorithm 680, collected algorithms from ACM. * Reference - GPM Poppe, CMJ Wijers: More efficient computation of * the complex error-function, ACM Trans. Math. Software. * Vol. 16, No. 1, Pp. 47. * * Modified by W. V. Snyder for inclusion in Math77: * Reorganize checking for overflow and loss of precision so * there are no redundant or unnecessary checks. In the process, * the range of applicability is expanded to the entire upper * half-plane. * Reorganize some calculations to be immune from overflow. * Split loop for two outer regions into two loops -- faster in * region. * Use D1MACH to fetch machine characteristics. * Use Math77 error message processor. * * Given a complex number z = (xi,yi), this subroutine computes * the value of the Faddeeva function w(z) = exp(-z**2)*erfc(-i*z), * where erfc is the complex complementary error function and i * means sqrt(-1). * The accuracy of the algorithm for z in the 1st and 2nd quadrant * is 14 significant digits; in the 3rd and 4th it is 13 significant * digits outside a circular region with radius 0.126 around a zero * of the function. * * * Argument list * Z [in] = real and imaginary parts of z in Z(1) and Z(2) * W [out] = real and imaginary parts of w(z) in W(1) and W(2) * FLAG [out] = an error flag indicating the status of the * computation. Type INTEGER, with values having the following * meaning: * 0 : No error condition, * -1 : Overflow would occur, * +1 : There would be no significant digits in the answer. * * * The routine is not underflow-protected but any variable can be * put to zero upon underflow. * *-- S version uses CWOFZ, R1MACH, serv1 *-- D version uses ZWOFZ, D1MACH, derv1 * * * RSQPD2 = 2/sqrt(pi) = reciprocal ( sqrt(pi) / 2 ). * LN2 = ln(2). * SQLGOV = sqrt(ln(RMAX)), where RMAX = the overflow limit for * floating point arithmetic. * LGOV2 = ln(RMAX) - ln(2) * LGUND = ln(RMIN), where RMIN = underflow limit. * MXGONI = the largest possible argument of sin or cos, restricted * here to sqrt ( pi / (2*round-off-limit) ). * INOGXM = 1 / MXGONI. * The reason these values are needed as defined * will be explained by comments in the code. * */ if (lgov2 <= 0.0e0) { lgov2 = log( DBL_MAX ); sqlgov = sqrt( lgov2 ); lgov2 -= LN2; lgund = log( DBL_MIN ); mxgoni = sqrt( 8.0e0/DBL_EPSILON )/RSQPD2; inogxm = 1.0e0/mxgoni; } xabs = fabs( Z[1] ); yabs = fabs( Z[2] ); x = xabs/6.3e0; y = yabs/4.4e0; if (x > y) { qrho = x*sqrt( 1.0e0 + powi(y/x,2) ); } else if (y == 0.0e0) { qrho = 0.0e0; } else { qrho = y*sqrt( 1.0e0 + powi(x/y,2) ); } a = qrho < 0.292e0; if (a) { /* qrho .lt. 0.292, equivalently qrho**2 .lt. 0.085264: the Fadeeva * function is evaluated using a power-series (Abramowitz and * Stegun, equation (7.1.5), p.297). * N is the minimum number of terms needed to obtain the required * accuracy. * * We know xquad and exp(-xquad) and yqyad and sin(yquad) won't * cause any trouble here, because qrho .lt. 1. */ xquad = (xabs - yabs)*(xabs + yabs); yquad = 2.0e0*xabs*yabs; n = (long)( 6.5e0 + 72.0e0*(1.0e0 - 0.85e0*y)*qrho ); j = 2*n + 1; xsum = RSQPD2/j; ysum = 0.0e0; for (i = n; i >= 1; i--) { j -= 2; w1 = (xsum*xquad - ysum*yquad)/i; ysum = (xsum*yquad + ysum*xquad)/i; xsum = w1 + RSQPD2/j; } u1 = 1.0e0 - (xsum*yabs + ysum*xabs); v1 = xsum*xabs - ysum*yabs; w1 = exp( -xquad ); u2 = w1*cos( yquad ); v2 = -w1*sin( yquad ); u = u1*u2 - v1*v2; v = u1*v2 + v1*u2; } else { rx = 0.0e0; ry = 0.0e0; sx = 0.0e0; sy = 0.0e0; /* The loops in both branches of the IF block below are similar. * They could be combined to reduce space, but extra tests and * unnecessary computation would be needed. * */ if (qrho < 1.0e0) { /* 0.292 .le. qrho .lt. 1.0: w(z) is evaluated by a truncated * Taylor expansion, where the Laplace continued fraction * is used to calculate the derivatives of w(z). * KAPN is the minimum number of terms in the Taylor expansion * needed to obtain the required accuracy. * NU is the minimum number of terms of the continued fraction * needed to calculate the derivatives with the required * accuracy. * x*x + y*y is more accurate than qrho*qrho here: */ c = (1.0e0 - y)*sqrt( 1.0e0 - x*x - y*y ); h = 1.88e0*c; h2 = 2.0e0*h; nu = (long)( 17.5e0 + 26.0e0*c ); kapn = (long)( 8.5e0 + 34.0e0*c ); /* Select kapn so qlamda doesn't underflow. Small kapn is good * (when possible) for performance also. */ if (h2 < 0.25e0) kapn = min( kapn, 1 + (long)( lgund/log( h2 ) ) ); qlamda = powi(h2,kapn - 1); /* 0 < qlamda < 3.76**41 < 3.85d23. */ for (n = nu; n >= (kapn + 1); n--) { tx = yabs + h + n*rx; ty = xabs - n*ry; /* No overflow because tx*rx + ty*ry = 1 and 0.292 < qrho < 1: */ c = 0.5e0/(tx*tx + ty*ty); rx = tx*c; ry = ty*c; } for (n = kapn; n >= 1; n--) { tx = yabs + h + n*rx; ty = xabs - n*ry; /* No overflow because tx*rx + ty*ry = 1 and 0.292 < qrho < 1: */ c = 0.5e0/(tx*tx + ty*ty); rx = tx*c; ry = ty*c; tx = qlamda + sx; sx = rx*tx - ry*sy; sy = ry*tx + rx*sy; qlamda /= h2; } u = RSQPD2*sx; v = RSQPD2*sy; } else { /* qrho .ge. 1.O: w(z) is evaluated using the Laplace continued * fraction. * NU is the minimum number of terms needed to obtain the * required accuracy. */ nu = (long)( 4.5e0 + (1442.0e0/(26.0e0*qrho + 77.0e0)) ); for (n = nu; n >= 1; n--) { tx = yabs + n*rx; ty = xabs - n*ry; if (tx > fabs( ty )) goto L_50; /* rx = 0.5*tx/(tx**2+ty**2) and ry = 0.5*ty/(tx**2+ty**2), * computed without overflow. Underflow is OK. */ c = tx/ty; ry = 0.5e0/(ty*(1.0e0 + c*c)); rx = ry*c; } goto L_60; /* Once tx>abs(ty), it stays that way. */ L_50: ; /* rx = 0.5*tx/(tx**2+ty**2) and ry = 0.5*ty/(tx**2+ty**2), * computed without overflow. Underflow is OK. */ c = ty/tx; rx = 0.5e0/(tx*(1.0e0 + c*c)); ry = rx*c; n -= 1; if (n == 0) goto L_60; tx = yabs + n*rx; ty = xabs - n*ry; goto L_50; L_60: u = RSQPD2*rx; v = RSQPD2*ry; } if (yabs == 0.0e0) { if (xabs > sqlgov) { u = 0.0e0; } else { u = exp( -SQ(xabs) ); } } } /* Evaluation of w(z) in the other quadrants. * */ if (Z[2] < 0.0e0) { if (a) { u2 += u2; v2 += v2; } else { /* Check whether sin(2*xabs*yabs) has any precision, without * allowing 2*xabs*yabs to overflow. */ if (yabs > xabs) { if (yabs > inogxm) { /* The following protects 2*exp(-z**2) against overflow. */ if (lgov2/yabs < yabs - xabs*(xabs/yabs)) goto L_100; w1 = 2.0e0*exp( (yabs - xabs)*(yabs + xabs) ); uv = fmin( fabs( u ), fabs( v ) ); if (w1 > uv) { if (xabs > mxgoni/yabs) goto L_110; } else { /* We put xabs*(w1/uv) here instead of simply xabs because * loss of precision in sin and cos will be diminished * relative to uv by w1. */ if (xabs*(w1/uv) > mxgoni/yabs) goto L_110; } } } else if (xabs > inogxm) { if (lgov2/xabs < xabs - yabs*(yabs/xabs)) { /* (yabs-xabs)*(yabs+xabs) might have overflowed, but in that * case, exp((yabs-xabs)*(yabs+xabs)) would underflow. */ u2 = 0.0e0; v2 = 0.0e0; goto L_80; } /* (yabs-xabs)*(yabs+xabs) can't overflow here. */ w1 = 2.0e0*exp( (yabs - xabs)*(yabs + xabs) ); uv = fmin( fabs( u ), fabs( v ) ); if (w1 > uv) { if (yabs > mxgoni/xabs) goto L_110; } else { /* We put yabs*(w1/uv) here instead of simply yabs because * loss of precision in sin and cos will be diminished * relative to uv by w1. */ if (yabs*(w1/uv) > mxgoni/xabs) goto L_110; } } yquad = 2.0e0*xabs*yabs; u2 = w1*cos( yquad ); v2 = -w1*sin( yquad ); L_80: ; } u = u2 - u; v = v2 - v; if (Z[1] > 0.0e0) v = -v; } else { if (Z[1] < 0.0e0) v = -v; } *flag = 0; W[1] = u; W[2] = v; return; /* Overflow */ L_100: *flag = -1; ermsg( "ZWOFZ", -1, 2, "EXP(-REAL(Z**2)) would overflow, Y .lt. - SQRT(X**2 + LN(infinity/2))" , ',' ); derv1( "LN(infinity/2)", lgov2, ',' ); goto L_120; /* No significant digits */ L_110: *flag = 1; ermsg( "ZWOFZ", 1, 2, "Too few significant digits in COS(IMAG(Z**2))," , ',' ); ermor( "2*Y*ABS(X) .lt. -SQRT(2*pi/eps)", ',' ); derv1( "SQRT(2*pi/eps)", mxgoni, ',' ); L_120: derv1( "REAL(Z)", Z[1], ',' ); derv1( "IMAG(Z)", Z[2], '.' ); W[1] = DBL_MAX; W[2] = W[1]; return; } /* end of function */