MINICERN/mathlib/gen/u/dtclgn64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:48  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       SUBROUTINE DTCLGN(J1,J2,J3,M1,M2,M3,DNUM,DDEN,KPEX)
  12       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  13 #endif
  14 #if !defined(CERNLIB_DOUBLE)
  15       SUBROUTINE DTCLGN(J1,J2,J3,M1,M2,M3,DNUM,DDEN,KPEX)
  16 #include "gen/imp64.inc"
  17 #endif
  18       CHARACTER NAME*(*)
  19       CHARACTER*80 ERRTXT
  20       PARAMETER (NAME = 'DTCLGN')
  21
  22       PARAMETER (MXP = 40)
  23
  24       DIMENSION KPEX(*)
  25       DIMENSION IPRIM(MXP),PRIM(MXP),LS(MXP,12),LR(MXP,100)
  26       DIMENSION ISD(11),JS(100),IA(7)
  27
  28       DATA (IA(I),I=1,7) /-1,1,1,-1,-1,1,1/
  29
  30       DATA (IPRIM(I),I=1,MXP)
  31      1/  2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
  32      2  59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,
  33      3 137,139,149,151,157,163,167,173/
  34
  35       DO 11 J = 1,MXP
  36       LS(J,12)=0
  37       PRIM(J)=IPRIM(J)
  38    11 KPEX(J)=0
  39       DNUM=0
  40       DDEN=1
  41       IF(J1 .LT. 0 .OR. J2 .LT. 0 .OR. J3 .LT. 0) GO TO 99
  42       IF(MOD(J1+J2+J3,2) .NE. 0) GO TO 99
  43       MX=ABS(M1)
  44       MY=ABS(M2)
  45       MZ=ABS(M3)
  46       IF(J1 .LT. MX .OR. J2 .LT. MY .OR. J3 .LT. MZ) GO TO 99
  47       IF(MOD(J1+MX,2) .NE. 0 .OR. MOD(J2+MY,2) .NE. 0) GOTO 99
  48       IF(MOD(J3+MZ,2) .NE. 0) GOTO 99
  49       IF(M1+M2 .NE. M3) GO TO 99
  50       ISD(1)=J3+1
  51       ISD(2)=(J3+J1-J2)/2
  52       ISD(3)=(J3-J1+J2)/2
  53       ISD(4)=(J1+J2-J3)/2
  54       ISD(5)=(J3-M3)/2
  55       ISD(6)=(J3+M3)/2
  56       ISD(7)=(J1+J2+J3)/2+1
  57       ISD(8)=(J1-M1)/2
  58       ISD(9)=(J1+M1)/2
  59       ISD(10)=(J2-M2)/2
  60       ISD(11)=(J2+M2)/2
  61
  62       IF(ISD(2) .LT. 0 .OR. ISD(3) .LT. 0 .OR. ISD(4) .LT. 0) GO TO 99
  63
  64       DO 41 K = 1,MXP
  65       DO 41 L = 1,11
  66    41 LS(K,L)=0
  67       IF(ISD(1) .GT. 1) THEN
  68        KOVF=1
  69        N1=ISD(1)
  70        K=0
  71    42  K=K+1
  72        IF(K .EQ. MXP+1) GO TO 48
  73    43  IF(MOD(N1,IPRIM(K)) .NE. 0) GO TO 42
  74        LS(K,1)=LS(K,1)+1
  75        N1=N1/IPRIM(K)
  76        IF(N1 .NE. 1) GO TO 43
  77    48  KOVF=MAX(KOVF,N1)
  78        IF(KOVF .NE. 1) THEN
  79         WRITE(ERRTXT,101) J1,J2,J3,M1,M2,M3
  80         CALL MTLPRT(NAME,'U112.1',ERRTXT)
  81         GO TO 99
  82        ENDIF
  83       ENDIF
  84       DO 59 L = 2,11
  85       DO 54 N = 2,ISD(L)
  86       KOVF=1
  87       N1=N
  88       K=0
  89    52 K=K+1
  90       IF(K .EQ. MXP+1) GO TO 58
  91    53 IF(MOD(N1,IPRIM(K)) .NE. 0) GO TO 52
  92       LS(K,L)=LS(K,L)+1
  93       N1=N1/IPRIM(K)
  94       IF(N1 .NE. 1) GO TO 53
  95    58 KOVF=MAX(KOVF,N1)
  96       IF(KOVF .NE. 1) THEN
  97        WRITE(ERRTXT,101) J1,J2,J3,M1,M2,M3
  98        CALL MTLPRT(NAME,'U112.1',ERRTXT)
  99        GO TO 99
 100       ENDIF
 101    54 CONTINUE
 102    59 CONTINUE
 103
 104       DO 5 J = 1,MXP
 105     5 LS(J,12)=LS(J,12)+LS(J,1)+LS(J,2)+LS(J,3)+LS(J,4)+LS(J,5)+LS(J,6)
 106      1        -LS(J,7)-LS(J,8)-LS(J,9)-LS(J,10)-LS(J,11)
 107       MM=0
 108       ISD(1)=J2+J3+M1
 109       ISD(2)=J1-M1
 110       ISD(3)=0
 111       ISD(4)=J3-J1+J2
 112       ISD(5)=J3+M3
 113       ISD(6)=J1-J2-M3
 114       ISD(7)=J2+M2
 115       DO 12 L = 1,7
 116    12 ISD(L)=ISD(L)/2-IA(L)
 117     3 DO 13 L = 1,7
 118    13 ISD(L)=ISD(L)+IA(L)
 119       IF(ISD(2) .LT. 0 .OR. ISD(6) .LT. 0) GO TO 3
 120       IF(ISD(1) .LT. 0 .OR. ISD(4) .LT. 0 .OR. ISD(5) .LT. 0) GO TO 4
 121
 122       DO 69 L = 1,6
 123       DO 61 K = 1,MXP
 124    61 LS(K,L)=0
 125       DO 64 N = 2,ISD(L)
 126       KOVF=1
 127       N1=N
 128       K=0
 129    62 K=K+1
 130       IF(K .EQ. MXP+1) GO TO 68
 131    63 IF(MOD(N1,IPRIM(K)) .NE. 0) GO TO 62
 132       LS(K,L)=LS(K,L)+1
 133       N1=N1/IPRIM(K)
 134       IF(N1 .NE. 1) GO TO 63
 135    68 KOVF=MAX(KOVF,N1)
 136       IF(KOVF .NE. 1) THEN
 137        WRITE(ERRTXT,101) J1,J2,J3,M1,M2,M3
 138        CALL MTLPRT(NAME,'U112.1',ERRTXT)
 139        GO TO 99
 140       ENDIF
 141    64 CONTINUE
 142    69 CONTINUE
 143
 144       MM=MM+1
 145       JS(MM)=(-1)**ISD(7)
 146       DO 9 J = 1,MXP
 147     9 LR(J,MM)=LS(J,1)+LS(J,2)-LS(J,3)-LS(J,4)-LS(J,5)-LS(J,6)
 148       GO TO 3
 149
 150     4 DO 31 J = 1,MXP
 151       LS(J,8)=10000
 152       DO 31 JK = 1,MM
 153    31 LS(J,8)=MIN(LS(J,8),LR(J,JK))
 154       BSUM=0
 155       DO 32 JM = 1,MM
 156       BNUM=1
 157       DO 33 J = 1,MXP
 158       JEX=LR(J,JM)-LS(J,8)
 159       IF(JEX .GT. 0) BNUM=BNUM*PRIM(J)**JEX
 160    33 CONTINUE
 161    32 BSUM=BSUM+BNUM*JS(JM)
 162       IF(BSUM .EQ. 0) GO TO 99
 163       DO 71 K = 1,MXP
 164    71 LS(K,1)=0
 165       ASUM=ABS(BSUM)
 166       IF(ASUM .GT. 1) THEN
 167        OVF=1
 168        A1=ASUM
 169        K=0
 170    72  K=K+1
 171        IF(K .EQ. MXP+1) GO TO 78
 172    73  IF(MOD(A1,PRIM(K)) .NE. 0) GO TO 72
 173        LS(K,1)=LS(K,1)+1
 174        A1=A1/PRIM(K)
 175        IF(A1 .NE. 1) GO TO 73
 176    78  OVF=MAX(OVF,A1)
 177        IF(OVF .NE. 1) THEN
 178         WRITE(ERRTXT,101) J1,J2,J3,M1,M2,M3
 179         CALL MTLPRT(NAME,'U112.1',ERRTXT)
 180         GO TO 99
 181        ENDIF
 182       ENDIF
 183       DO 22 J = 1,MXP
 184    22 KPEX(J)=LS(J,12)+2*(LS(J,8)+LS(J,1))
 185       BNUM=1
 186       BDEN=1
 187       DO 16 J = 1,MXP
 188       PP=PRIM(J)**ABS(KPEX(J))
 189       IF(KPEX(J) .GE. 0) THEN
 190        BNUM=BNUM*PP
 191       ELSE
 192        BDEN=BDEN*PP
 193       ENDIF
 194    16 CONTINUE
 195       DNUM=SIGN(BNUM,BSUM)
 196       DDEN=BDEN
 197    99 RETURN
 198   101 FORMAT('LIST OF PRIME NUMBERS EXHAUSTED',2X,3I5,2X,3I5)
 199       END