github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zsqrt.f (about)

     1        SUBROUTINE ZSQRT(AR, AI, BR, BI)
     2  C***BEGIN PROLOGUE  ZSQRT
     3  C***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
     4  C
     5  C     DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A)
     6  C
     7  C***ROUTINES CALLED  ZABS
     8  C***END PROLOGUE  ZSQRT
     9        DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DRT
    10        DOUBLE PRECISION ZABS
    11  
    12        DATA DRT , DPI / 7.071067811865475244008443621D-1,
    13       1                 3.141592653589793238462643383D+0/
    14        ZM = ZABS(CMPLX(AR,AI,kind=KIND(1.0D0)))
    15        ZM = DSQRT(ZM)
    16        IF (AR.EQ.0.0D+0) GO TO 10
    17        IF (AI.EQ.0.0D+0) GO TO 20
    18        DTHETA = DATAN(AI/AR)
    19        IF (DTHETA.LE.0.0D+0) GO TO 40
    20        IF (AR.LT.0.0D+0) DTHETA = DTHETA - DPI
    21        GO TO 50
    22     10 IF (AI.GT.0.0D+0) GO TO 60
    23        IF (AI.LT.0.0D+0) GO TO 70
    24        BR = 0.0D+0
    25        BI = 0.0D+0
    26        RETURN
    27     20 IF (AR.GT.0.0D+0) GO TO 30
    28        BR = 0.0D+0
    29        BI = DSQRT(DABS(AR))
    30        RETURN
    31     30 BR = DSQRT(AR)
    32        BI = 0.0D+0
    33        RETURN
    34     40 IF (AR.LT.0.0D+0) DTHETA = DTHETA + DPI
    35     50 DTHETA = DTHETA*0.5D+0
    36        BR = ZM*DCOS(DTHETA)
    37        BI = ZM*DSIN(DTHETA)
    38        RETURN
    39     60 BR = ZM*DRT
    40        BI = ZM*DRT
    41        RETURN
    42     70 BR = ZM*DRT
    43        BI = -ZM*DRT
    44        RETURN
    45        END