1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2c Copyright (C) Bruno Pincon 3c 4c Copyright (C) 2012 - 2016 - Scilab Enterprises 5c 6c This file is hereby licensed under the terms of the GNU GPL v2.0, 7c pursuant to article 5.3.4 of the CeCILL v.2.1. 8c This file was originally licensed under the terms of the CeCILL v2.1, 9c and continues to be available under such terms. 10c For more information, see the COPYING file which you should have received 11c along with this program. 12 13 subroutine wacos(zr, zi, ar, ai) 14* 15* PURPOSE 16* Compute the arccosine of a complex number 17* a = ar + i ai = acos(z), z = zr + i zi 18* 19* CALLING LIST / PARAMETERS 20* subroutine wacos(zr,zi,ar,ai) 21* double precision zr,zi,ar,ai 22* 23* zr,zi: real and imaginary parts of the complex number 24* ar,ai: real and imaginary parts of the result 25* ar,ai may have the same memory cases than zr et zi 26* 27* REFERENCE 28* This is a Fortran-77 translation of an algorithm by 29* T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which 30* appears in their article : 31* "Implementing the Complex Arcsine and Arccosine 32* Functions Using Exception Handling", ACM, TOMS, 33* Vol 23, No. 3, Sept 1997, p. 299-335 34* 35* with some modifications so as don't rely on ieee handle 36* trap functions. 37* 38* AUTHOR 39* Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr> 40* Thanks to Tom Fairgrieve 41* 42 implicit none 43 44* PARAMETERS 45 double precision zr, zi, ar, ai 46 47* EXTERNAL FUNCTIONS 48 double precision dlamch, logp1 49 external dlamch, logp1 50 integer isanan 51 external isanan 52 53* CONSTANTS 54 double precision LN2, PI, HALFPI, Across, Bcross 55 parameter (LN2 = 0.6931471805599453094172321d0, 56 $ HALFPI = 1.5707963267948966192313216d0, 57 $ PI = 3.1415926535897932384626433d0, 58 $ Across = 1.5d0, 59 $ Bcross = 0.6417d0) 60* LOCAL VARIABLES 61 double precision x, y, A, B, R, S, Am1, szr, szi 62 63 64* STATIC VARIABLES 65 double precision LSUP, LINF, EPSM 66 save LSUP, LINF, EPSM 67 logical first 68 save first 69 data first /.true./ 70 71* TEXT 72* got f.p. parameters used by the algorithm 73 if (first) then 74 LSUP = sqrt(dlamch('o'))/8.d0 75 LINF = 4.d0*sqrt(dlamch('u')) 76 EPSM = sqrt(dlamch('e')) 77 first = .false. 78 endif 79 80* avoid memory pb ... 81 x = abs(zr) 82 y = abs(zi) 83 szr = sign(1.d0,zr) 84 szi = sign(1.d0,zi) 85 86 87 if (LINF .le. min(x,y) .and. max(x,y) .le. LSUP ) then 88* we are in the safe region 89 R = sqrt((x+1.d0)**2 + y**2) 90 S = sqrt((x-1.d0)**2 + y**2) 91 A = 0.5d0*(R + S) 92 B = x/A 93 94* compute the real part 95 if ( B .le. Bcross ) then 96 ar = acos(B) 97 elseif ( x .le. 1.d0 ) then 98 ar = atan( sqrt( 0.5d0*(A+x) * 99 $ ( (y**2)/(R+(x+1.d0)) + (S+(1.d0-x)) ) ) / x ) 100 else 101 ar = atan((y*sqrt(0.5d0*((A+x)/(R+(x+1.d0)) 102 $ +(A+x)/(S+(x-1.d0))))) / x) 103 endif 104 105* compute the imaginary part 106 if ( A .le. Across ) then 107 if ( x .lt. 1.d0 ) then 108 Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) 109 else 110 Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) 111 endif 112 ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) 113 else 114 ai = log(A + sqrt(A**2 - 1.d0)) 115 endif 116 117 else 118* HANDLE BLOC : evaluation in the special regions ... 119 if ( y .le. EPSM*abs(x-1.d0) ) then 120 if (x .lt. 1.d0 ) then 121 ar = acos(x) 122 ai = y/sqrt((1.d0+x)*(1.d0-x)) 123 else 124 ar = 0.d0 125 if ( x .le. LSUP ) then 126 ai = logp1((x-1.d0) + sqrt((x-1.d0)*(x+1.d0))) 127 else 128 ai = LN2 + log(x) 129 endif 130 endif 131 elseif (y .lt. LINF) then 132 if (isanan(x).eq.1) then 133 ar = x 134 ai = y 135 else 136 ar = sqrt(y) 137 ai = ar 138 endif 139 elseif (EPSM*y - 1.d0 .ge. x) then 140 ar = HALFPI 141 ai = LN2 + log(y) 142 elseif (x .gt. 1.d0) then 143 ar = atan(y/x) 144 ai = LN2 + log(y) + 0.5d0*logp1((x/y)**2) 145 else 146 if (isanan(x).eq.1) then 147 ar = x 148 else 149 ar = HALFPI 150 endif 151 A = sqrt(1.d0 + y**2) 152 ai = 0.5d0*logp1(2.d0*y*(y+A)) 153 endif 154 endif 155 156* recover the signs 157 if (szr .lt. 0.d0) then 158 ar = PI - ar 159 endif 160 161 if (y.ne.0.d0 .or. szr.lt.0.d0) then 162 ai = - szi * ai 163 endif 164 165 end 166 167 168