MINICERN/mathlib/gen/divon/quasi.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:03:24  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE QUASI (X,NDIM,NGIVE,NF)
  11 C--      CORRECTED 830401 CERN GEN PAM 1.05
  12 C--      INITIAL VALUE OF I CHANGED FROM 1 TO 2 FOR CASE OF NPT LT 7
  13 C--      CORRECTED 850503 CERN GEN PAM 1.16
  14 C--      17'TH IN LIST OF PRIME NUMBERS CORRECTED FROM 1197 TO 1187
  15       INTEGER NDIM,N,NGIVE
  16       REAL X(10,NGIVE)
  17       INTEGER PRIME(21),AA(20,9)
  18       REAL THETA(10)
  19       INTEGER D0,N0,A,P,IPT
  20       DATA D0,N0 /2*0/
  21       DATA PRIME / 7,19,29,37,47,97,149,199,293,397,499,599,691,797,887,
  22      1997,1187, 1499,1789,1999,9999999 /
  23       DATA AA / 2,7,12,27,29,61,44,76,81,163,209,165,390,346,192,705, 65
  24      16,629,1037,878, 2,3,20,13,30,78,67,73,103,274,405,248,585,305,674,
  25      2650, 358,526,1540,1082, 2,9,21,16,37,64,29,114,211,157,241,472,371
  26      3,477,113,252, 736,1205,853,1486, 2,2,27,35,44,8,110,42,258,104,445
  27      4,21,640,665,678,535, 1003,56,1634,1781, 2,17,15,29,38,76,85,183,28
  28      5,356,452,522,655,787,841,697, 1073,1025,1025,962, 2,17,2,2,2,67,12
  29      68,94,257,393,395,395,355,341,397,361, 383,367,383,398, 2,17,2,19,2
  30      7,79,139,149,170,165,169,161,157,171,149,167, 157,171,133,151, 2,17
  31      8,2,19,45,2,139,67,170,87,169,85,157,89,81,167,157,75,89,83, 2,17,2
  32      97,35,45,47,51,45,170,87,53,51,157,89,81,167,157,53,53,53 /
  33       N=ABS(NF)
  34       IF(NGIVE.LT.1.OR.NGIVE.GT.N) NGIVE=N
  35       IF(N.EQ.N0 .AND. NDIM.EQ.D0) GOTO 60
  36       N0=N
  37       D0=NDIM
  38       IPT=0
  39 C--      INITIAL VALUE CHANGED FROM 1 TO 2, PAM 1.05
  40       I=2
  41       GOTO 20
  42  10   I=I+1
  43  20   IF((I).GT.(21)) GOTO 30
  44       IF(PRIME(I).LE.N) GOTO 10
  45  30   P=PRIME(I-1)
  46       A=AA(I-1,NDIM-1)
  47       THETA(1)=1.0E+0
  48       THETA(2)=A
  49       DO 40 I=3,NDIM
  50       THETA(I)=MOD(THETA(2)*THETA(I-1),REAL(P))
  51  40   CONTINUE
  52       DO 50 I=1,NDIM
  53       THETA(I)=THETA(I)/P
  54  50   CONTINUE
  55  60   IF(IPT+NGIVE.GT.P) NGIVE=P-IPT
  56       IF(NF.LT.0) RETURN
  57       DO 80 K=1,NGIVE
  58       IPT=IPT+1
  59       DO 70 I=1,NDIM
  60       TERM=IPT*THETA(I)
  61       X(I,K)=ABS(2.0E+0*MOD(TERM,1.0E+0)-1.0E+0)
  62  70   CONTINUE
  63  80   CONTINUE
  64       IF(IPT.EQ.P) IPT=0
  65       RETURN
  66       END