MINICERN/mathlib/gen/c/rteq464.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:53  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       SUBROUTINE RRTEQ4(A,B,C,D,Z,DC,MT)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       SUBROUTINE DRTEQ4(A,B,C,D,Z,DC,MT)
  15 #endif
  16 #include "gen/imp64.inc"
  17 #include "gen/defc64.inc"
  18      +  I,Z,Z0,W1,W2,W3
  19 #include "gen/def128.inc"
  20      +      ZQ1,AA,PP,QQ,RR,Q1,Q2,Q3,Q4,Q8
  21       DIMENSION Z(*),Z0(4),U(3),V(3)
  22
  23       PARAMETER(I = (0,1), ZD1 = 1, ZQ1 = 1)
  24       PARAMETER(R4 = ZD1/4, R12 = ZD1/12)
  25       PARAMETER(Q2 = ZQ1/2, Q4 = ZQ1/4, Q8 = ZQ1/8)
  26       PARAMETER(Q1 = 3*ZQ1/8, Q3 = 3*ZQ1/16)
  27
  28       IF(B .EQ. 0 .AND. C .EQ. 0) THEN
  29        IF(D .EQ. 0) THEN
  30         MT=1
  31         Z(1)=-A
  32         Z(2)=0
  33         Z(3)=0
  34         Z(4)=0
  35         DC=0
  36         RETURN
  37        ELSEIF(A .EQ. 0) THEN
  38         IF(D .GT. 0) THEN
  39          MT=2
  40          Z(1)=SQRT(I*SQRT(D))
  41          Z(2)=-Z(1)
  42          Z(4)=SQRT(-Z(1)**2)
  43          Z(3)=-Z(4)
  44         ELSE
  45          MT=3
  46          Z(1)=SQRT(SQRT(-D))
  47          Z(2)=-Z(1)
  48          Z(3)=SQRT(-Z(1)**2)
  49          Z(4)=-Z(3)
  50         ENDIF
  51         DC=(-R12*D)**3
  52         RETURN
  53        ENDIF
  54       ENDIF
  55       AA=A**2
  56       PP=B-Q1*AA
  57       QQ=C-Q2*A*(B-Q4*AA)
  58       RR=D-Q4*(A*C-Q4*AA*(B-Q3*AA))
  59       RC=Q2*PP
  60       SC=Q4*(Q4*PP**2-RR)
  61       TC=-(Q8*QQ)**2
  62 #if defined(CERNLIB_DOUBLE)
  63       CALL DRTEQ3(RC,SC,TC,U,DC)
  64 #endif
  65 #if !defined(CERNLIB_DOUBLE)
  66       CALL RRTEQ3(RC,SC,TC,U,DC)
  67 #endif
  68       Q=QQ
  69       H=R4*A
  70       IF(DC .EQ. 0) U(3)=U(2)
  71       IF(DC .LE. 0) THEN
  72        MT=2
  73        V(1)=ABS(U(1))
  74        V(2)=ABS(U(2))
  75        V(3)=ABS(U(3))
  76        V1=MAX(V(1),V(2),V(3))
  77        IF(V1 .EQ. V(1)) THEN
  78         K1=1
  79         V2=MAX(V(2),V(3))
  80        ELSEIF(V1 .EQ. V(2)) THEN
  81         K1=2
  82         V2=MAX(V(1),V(3))
  83        ELSE
  84         K1=3
  85         V2=MAX(V(1),V(2))
  86        ENDIF
  87        IF(V2 .EQ. V(1)) THEN
  88         K2=1
  89        ELSEIF(V2 .EQ. V(2)) THEN
  90         K2=2
  91        ELSE
  92         K2=3
  93        ENDIF
  94        W1=SQRT(U(K1)+I*0)
  95        W2=SQRT(U(K2)+I*0)
  96       ELSE
  97        MT=3
  98        W1=SQRT(U(2)+I*U(3))
  99        W2=SQRT(U(2)-I*U(3))
 100       ENDIF
 101       W3=0
 102       IF(W1*W2 .NE. 0) W3=-Q/(8*W1*W2)
 103       Z0(1)=W1+W2+W3-H
 104       Z0(2)=-W1-W2+W3-H
 105       Z0(3)=-W1+W2-W3-H
 106       Z0(4)=W1-W2-W3-H
 107       IF(MT .EQ. 2) THEN
 108        IF(U(K1) .GE. 0 .AND. U(K2) .GE. 0) THEN
 109         MT=1
 110         DO 1 J = 1,4
 111         RZ0=Z0(J)
 112     1   Z(J)=RZ0
 113        ELSEIF(U(K1) .GE. 0 .AND. U(K2) .LT. 0) THEN
 114         Z(1)=Z0(1)
 115         Z(2)=Z0(4)
 116         Z(3)=Z0(3)
 117         Z(4)=Z0(2)
 118        ELSEIF(U(K1) .LT. 0 .AND. U(K2) .GE. 0) THEN
 119         Z(1)=Z0(1)
 120         Z(2)=Z0(3)
 121         Z(3)=Z0(4)
 122         Z(4)=Z0(2)
 123        ELSEIF(U(K1) .LT. 0 .AND. U(K2) .LT. 0) THEN
 124         Z(1)=Z0(1)
 125         Z(2)=Z0(2)
 126         Z(3)=Z0(4)
 127         Z(4)=Z0(3)
 128        ENDIF
 129       ELSEIF(MT .EQ. 3) THEN
 130        DO 2 J = 1,2
 131        RZ0=Z0(J)
 132     2  Z(J)=RZ0
 133        Z(3)=Z0(4)
 134        Z(4)=Z0(3)
 135       ENDIF
 136       RETURN
 137       END