c Routine AS 93 returns frequencies.   The following short routine calculates
c the distribution function from these frequencies (overwriting them).
c The calling arguments are as for AS 93.   The distribution function is
c returned in array A1.   The first element in A1 is F(ASTART).   N.B. ASTART
c is a real variable.
c
      subroutine wprob(test, other, astart, a1, l1, a2, a3, ifault)
      integer test, other, l1, ifault
      real astart, a1(l1), a2(l1), a3(l1)
c
c     Local variables
c
      real zero, sum
      data zero /0.0/
c
      call gscale(test, other, astart, a1, l1, a2, a3, ifault)
      if (ifault .ne. 0) return
c
c     Scale column of F
c
      nrows = 1 + (test * other)/2
      sum = zero
      do 10 i = 1, nrows
	sum = sum + a1(i)
	a1(i) = sum
   10 continue
      do 20 i = 1, nrows
   20 a1(i) = a1(i) / sum
c
      return
      end

c----------------------------------------------------------------------

	SUBROUTINE GSCALE(TEST, OTHER, ASTART, A1, L1, A2, A3, IFAULT)
C
C	  ALGORITHM AS 93 APPL. STATIST. (1976) VOL.25, NO.1
C
C	  FROM THE SIZES OF TWO SAMPLES THE DISTRIBUTION OF THE
C	  ANSARI-BRADLEY TEST FOR SCALE IS GENERATED IN ARRAY A1.
C
	REAL ASTART, A1(L1), A2(L1), A3(L1), AI, ONE, FPOINT
	INTEGER TEST, OTHER
	LOGICAL SYMM
	DATA ONE /1.0/
C
C	  TYPE CONVERSION (EFFECT DEPENDS ON TYPE STATEMENT ABOVE).
C
	FPOINT(I) = I
C
C	  CHECK PROBLEM SIZE AND DEFINE BASE VALUE OF THE DISTRIBUTION.
C
	M = MIN0(TEST, OTHER)
	IFAULT = 2
	IF (M. LT. 0) RETURN
	ASTART = FPOINT((TEST + 1) / 2) * FPOINT(1 + TEST / 2)
	N = MAX0(TEST, OTHER)
C
C	  CHECK SIZE OF RESULT ARRAY.
C
	IFAULT = 1
	LRES = 1 + (M * N) / 2
	IF (L1 .LT. LRES) RETURN
	SYMM = MOD(M + N, 2) .EQ. 0
C
C	  TREAT SMALL SAMPLES SEPARATELY.
C
	MM1 = M - 1
	IF (M .GT. 2) GOTO 5
C
C	  START-UP PROCEDURES ONLY NEEDED.
C
	IF (MM1) 1, 2, 3
C
C	  ONE SAMPLE ONLY.
C
1	A1(1) = ONE
	GOTO 15
C
C	  SMALLER SAMPLE SIZE = 1.
C
2	 CALL START1(N, A1, L1, LN1)
	GOTO 4
C
C	  SMALLER SAMPLE SIZE = 2.
C
3	CALL START2(N, A1, L1, LN1)
C
C	RETURN IF A1 IS NOT IN REVERSE ORDER.
C
4	IF (SYMM .OR. (OTHER .GT. TEST)) GOTO 15
	GOTO 13
C
C	  FULL GENERATOR NEEDED
C	  SET UP INITIAL CONDITIONS (DEPENDS ON MOD(N, 2)).
C
5	NM1 = N - 1
	NM2 = N - 2
	MNOW = 3
	NC = 3
	IF (MOD(N, 2) .EQ. 1) GOTO 6
C	  SET UP FOR EVEN N.
C
	N2B1 = 3
	N2B2 = 2
	CALL START2(N, A1, L1, LN1)
	CALL START2(NM2, A3, L1, LN3)
	CALL START1(NM1, A2, L1, LN2)
	GOTO 8
C
C	  SET UP FOR ODD N.
C
6	N2B1 = 2
	N2B2 = 3
	CALL START1(N, A1, L1, LN1)
	CALL START2(NM1, A2, L1, LN2)
C
C	  INCREASE ORDER OF DISTRIBUTION IN A1 BY 2
C	  (USING A2 AND IMPLYING A3).
C
7	CALL FRQADD(A1, LN1, L1OUT, L1, A2, LN2, N2B1)
	LN1 = LN1 + N
	CALL IMPLY(A1, L1OUT, LN1, A3, LN3, L1, NC)
	NC = NC + 1
	IF (MNOW .EQ. M) GOTO 9
	MNOW = MNOW + 1
C
C	  INCREASE ORDER OF DISTRIBUTION IN A2 BY 2 (USING A3).
C
8	CALL FRQADD(A2, LN2, L2OUT, L1, A3, LN3, N2B2)
	LN2 = LN2 + NM1
	CALL IMPLY(A2, L2OUT, LN2, A3, J, L1, NC)
	NC = NC + 1
	IF (MNOW .EQ. M) GOTO 9
	MNOW = MNOW + 1
	GOTO 7
C
C	  IF SYMMETRICAL, RESULTS IN A1 ARE COMPLETE.
C
9	IF (SYMM) GOTO 15
C
C	  FOR A SKEW RESULT ADD A2 (OFFSET) INTO A1.
C
	KS = (M + 3) / 2
	J = 1
	DO 12 I = KS, LRES
	IF (I .GT. LN1) GOTO 10
	A1(I) = A1(I) + A2(J)
	GOTO 11
10	A1(I) = A2(J)
11	J = J + 1
12	CONTINUE
C
C	  DISTRIBUTION IN A1 POSSIBLY IN REVERSE ORDER.
C
	IF (OTHER .LT. TEST) GOTO 15
C
C	  REVERSE THE RESULTS IN A1.
C
13	J = LRES
	NDO = LRES / 2
	DO 14 I = 1, NDO
	AI = A1(I)
	A1(I) =A1(J)
	A1(J) = AI
	J = J - 1
14	CONTINUE
C
C 	  FINAL RESULTS NOW IN A1.
C
15	IFAULT = 0
	RETURN
	END

	SUBROUTINE START1(N, F, L, LOUT)
C
C	  ALGORITHM AS 93.1 APPL. STATIST. (1976) VOL.25, NO.1
C
C	  GENERATES A 1,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
	REAL F(L), ONE, TWO
	DATA ONE, TWO /1.0, 2.0/
	LOUT = 1 + N / 2
	DO 1 I = 1, LOUT
1	F(I) = TWO
	IF (MOD(N, 2) .EQ. 0) F(LOUT) = ONE
	RETURN
	END
C
	SUBROUTINE START2(N, F, L, LOUT)
C
C	  ALGORITHM AS 93.2 APPL. STATIST. (1976) VOL.25, NO.1
C
C	  GENERATES A 2,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
	REAL F(L), ONE, TWO, THREE, FOUR
	DATA ONE, TWO, THREE, FOUR /1.0, 2.0, 3.0, 4.0/
C
C	  DERIVE F FOR 2, NU, WHERE NU IS HIGHEST EVEN INTEGER
C	  LESS THAN OR EQUAL TO N.
C	  DEFINE NU AND ARRAY LIMITS.
C
	NU = N - MOD(N, 2)
	J = NU + 1
	LOUT = J
	LT1 = LOUT + 1
	NDO = LT1 / 2
	A = ONE
	B = THREE
C
C	  GENERATE THE SYMMETRICAL 2,NU DISTRIBUTION.
C
	DO 1 I = 1, NDO
	F(I) = A
	F(J) = A
	J = J - 1
	A = A + B
	B = FOUR - B
1	CONTINUE
	IF (NU .EQ. N) RETURN
C
C	  ADD AN OFFSET 1,N DISTRIBUTION INTO F TO GIVE 2,N RESULT.
C
	NU = NDO + 1
	DO 2 I = NU, LOUT
2	F(I) = F(I) + TWO
	F(LT1) = TWO
	LOUT = LT1
	RETURN
	END
C
	SUBROUTINE FRQADD(F1, L1IN, L1OUT, L1, F2, L2, NSTART)
C
C	  ALGORITHM AS 93.3 APPL. STATIST. (1976) VOL.25, NO.1
C
C	  ARRAY F1 HAS TWICE THE CONTENTS OF ARRAY F2 ADDED INTO IT
C	  STARTING WITH ELEMENTS NSTART AND 1 IN F1 AND F2 RESPECTIVELY.
C
	REAL F1(L1), F2(L2), MUL2
	DATA MUL2 /2.0/
	I2 = 1
	DO 1 I1 = NSTART, L1IN
	F1(I1) = F1(I1) + MUL2 * F2(I2)
	I2 = I2 + 1
1	CONTINUE
	NXT = L1IN + 1
	L1OUT = L2 + NSTART - 1
	DO 2 I1 = NXT, L1OUT
	F1(I1) = MUL2 * F2(I2)
	I2 = I2 + 1
2	CONTINUE
	NSTART = NSTART + 1
	RETURN
	END
C
	SUBROUTINE IMPLY(F1, L1IN, L1OUT, F2, L2, L2MAX, NOFF)
C
C	  ALGORITHM AS 93.4 APPL. STATIST. (1976) VOL.25, NO.1
C
C	  GIVEN L1IN ELEMENTS OF AN ARRAY F1, A SYMMETRICAL
C	  ARRAY F2 IS DERIVED AND ADDED ONTO F1, LEAVING THE
C	  FIRST NOFF ELEMENTS OF F1 UNCHANGED AND GIVING A
C	  SYMMETRICAL RESULT OF L1OUT ELEMENTS IN F1.
C
	REAL F1(L1OUT), F2(L2MAX), SUM, DIFF
C
C	  SET-UP SUBSCRIPTS AND LOOP COUNTER.
C
	I2 = 1 - NOFF
	J1 = L1OUT
	J2 = L1OUT - NOFF
	L2 = J2
	J2MIN = (J2 + 1) / 2
	NDO = (L1OUT + 1) / 2
C
C	  DERIVE AND IMPLY NEW VALUES FROM OUTSIDE INWARDS.
C
	DO 6 I1 = 1, NDO
C
C	  GET NEW F1 VALUE FROM SUM OF L/H ELEMENTS OF
C	  F1 + F2 (IF F2 IS IN RANGE).
C
	IF (I2 .GT. 0) GOTO 1
	SUM = F1(I1)
	GOTO 2
1	SUM = F1(I1) + F2(I2)
C
C	  REVISE LEFT ELEMENT OF F1.
C
	F1(I1) = SUM
C
C	  IF F2 NOT COMPLETE IMPLY AND ASSIGN F2 VALUES
C	  AND REVISE SUBSCRIPTS.
C
2	I2 = I2 + 1
	IF (J2 .LT. J2MIN) GOTO 5
	IF (J1 .LE. L1IN) GOTO 3
	DIFF = SUM
	GOTO 4
3	DIFF = SUM - F1(J1)
4	F2(I1) = DIFF
	F2(J2) = DIFF
	J2 = J2 - 1
C
C	  ASSIGN R/H ELEMENT OF F1 AND REVISE SUBSCRIPT.
C
5	F1(J1) = SUM
	J1 = J1 - 1
6	CONTINUE
	RETURN
	END