MINICERN/mathlib/gen/c/rmullz64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:52  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       SUBROUTINE DMULLZ(A,N,MAXITR,Z0)
  12 #endif
  13 #if !defined(CERNLIB_DOUBLE)
  14       SUBROUTINE RMULLZ(A,N,MAXITR,Z0)
  15 #endif
  16 #include "gen/imp64.inc"
  17 #include "gen/defc64.inc"
  18      +  Z0,Z,DX,X,X3,Y1,Y2,Y,TE(7)
  19       LOGICAL LSW
  20       DIMENSION A(0:*),Z0(*)
  21       CHARACTER*(*) NAME
  22 #if defined(CERNLIB_DOUBLE)
  23       PARAMETER( NAME='DMULLZ')
  24 #endif
  25 #if !defined(CERNLIB_DOUBLE)
  26       PARAMETER( NAME='RMULLZ')
  27 #endif
  28       PARAMETER (Z1 = 1, HF = Z1/2)
  29       PARAMETER (ETA1 = 1D-14, ETA2 = 6D-8, BIG = 1D20)
  30       PARAMETER (C1 = 0.9D0, C2 = 0.002D0, C3 = 0.1D0)
  31
  32 #include "gen/gcmpfun.inc"
  33       SUMABS(Z)=ABS(GREAL(Z))+ABS(GIMAG(Z))
  34
  35       IF(N .LE. 0) RETURN
  36       IF(A(0) .EQ. 0) THEN
  37        CALL MTLPRT(NAME,'C202.1','A(0) = 0')
  38        RETURN
  39       ENDIF
  40
  41       N1=N
  42     2 IF(N1 .EQ. 1) THEN
  43        Z0(1)=-A(1)/A(0)
  44        RETURN
  45       ENDIF
  46       IF(A(N1) .EQ. 0) THEN
  47        Z0(N1)=0
  48        N1=N1-1
  49        GO TO 2
  50       ENDIF
  51
  52       SCALE=ABS(A(N1)/A(0))**(Z1/N1)
  53       B=A(0)
  54       DO 6 I = 1,N1
  55       B=B*SCALE
  56     6 Z0(I)=A(I)/B
  57       IF(N1 .EQ. 2) GO TO 1
  58
  59    10 LSW=.TRUE.
  60       Y1=Z0(1)+1
  61       Y2=Z0(1)-1
  62       DO 11 I = 2,N1
  63       Y1=Z0(I)+Y1
  64    11 Y2=Z0(I)-Y2
  65       Y=Z0(N1)
  66       X=0
  67       DX=1
  68
  69       TE(1)=-2
  70    12 TE(2)=Y2/Y
  71       TE(3)=(Y1-Y2)/(Y*TE(1))
  72       DO 17 ITER = 1,MAXITR
  73       TE(4)=TE(2)-1
  74       TE(5)=(TE(4)-TE(3))/(TE(1)+1)
  75       TE(6)=HF*(TE(5)+TE(4))
  76       TE(7)=SQRT(TE(6)**2+TE(5))
  77       TE(1)=TE(6)+TE(7)
  78       TE(7)=TE(6)-TE(7)
  79       IF(ABS(TE(1)) .LE. ABS(TE(7))) THEN
  80        IF(TE(7) .EQ. 0) TE(7)=C1
  81        TE(1)=TE(7)
  82       ENDIF
  83       DX=DX/TE(1)
  84       X=DX+X
  85       EPSI=SUMABS(X)*ETA1
  86       IF(SUMABS(DX) .LT. EPSI .AND. SUMABS(Y) .LT. C2) GO TO 18
  87       Y2=Y
  88       Y=X+Z0(1)
  89       DO 21 I = 2,N1
  90    21 Y=Y*X+Z0(I)
  91       IF(Y .EQ. 0) GO TO 18
  92
  93    15 IF(SUMABS(Y) .GE. 100*SUMABS(Y2) .AND. SUMABS(DX) .GE. EPSI) THEN
  94        TE(1)=2*TE(1)
  95        DX=HF*DX
  96        X=X-DX
  97        Y=X+Z0(1)
  98        DO 22 I = 2,N1
  99    22  Y=Y*X+Z0(I)
 100        IF(Y .EQ. 0) GO TO 18
 101        GOTO 15
 102       ENDIF
 103
 104       TE(2)=Y2/Y
 105       TE(3)=(TE(2)/TE(1))*TE(4)
 106    17 CONTINUE
 107
 108       CN=ABS(Z0(N1))
 109       IF(ABS(CN-1) .LT. C3) THEN
 110        DO 40 I = 1,N1
 111    40  Z0(I)=BIG
 112        CALL MTLPRT(NAME,'C202.2','TOO MANY ITERATIONS')
 113        RETURN
 114       ENDIF
 115       S=CN**(-Z1/N1)
 116       SCALE=SCALE/S
 117       B=1
 118       DO 30 I = 1,N1
 119       B=B*S
 120    30 Z0(I)=Z0(I)*B
 121       GO TO 10
 122
 123    20 IF(ABS(GIMAG(X)) .LT. ABS(GREAL(X))*ETA2) GO TO 10
 124       LSW=.FALSE.
 125       X3=GCONJG(X)
 126       DX=GCONJG(DX)
 127       TE(1)=GCONJG(TE(1))
 128       X=X3-DX
 129       Y=X+Z0(1)
 130       DO 51 I = 2,N1
 131    51 Y=Y*X+Z0(I)
 132       IF(Y .EQ. 0) GO TO 18
 133       Y2=Y
 134       X=X-DX*TE(1)
 135       Y=X+Z0(1)
 136       DO 52 I = 2,N1
 137    52 Y=Y*X+Z0(I)
 138       IF(Y .EQ. 0) GO TO 18
 139       Y1=Y
 140       X=X3
 141       Y=X+Z0(1)
 142       DO 53 I = 2,N1
 143    53 Y=Y*X+Z0(I)
 144       IF(Y .NE. 0) GO TO 12
 145
 146    18 Z0(N1)=X*SCALE
 147       N1=N1-1
 148       Z0(1)=X+Z0(1)
 149       DO 19 I = 2,N1
 150    19 Z0(I)=Z0(I-1)*X+Z0(I)
 151       IF(N1 .GT. 2) THEN
 152        IF(LSW) GO TO 20
 153        GO TO 10
 154       ENDIF
 155
 156     1 Z0(2)=(SQRT((HF*Z0(1))**2-Z0(2))-HF*Z0(1))*SCALE
 157       Z0(1)=-Z0(1)*SCALE-Z0(2)
 158       RETURN
 159       END