MODULE Nonc_chisq
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL (dp), SAVE  :: aintl, ersm
REAL, SAVE       :: pi, aln28, sigsq, almax, almin, amean, c
INTEGER, SAVE    :: icount, ir, lim
LOGICAL, SAVE    :: ndtsrt, fail


CONTAINS


FUNCTION iround(x) RESULT(ival)
IMPLICIT NONE

REAL, INTENT(IN)  :: x
INTEGER           :: ival

REAL, PARAMETER  :: half = 0.5

ival = INT(x + SIGN(half,x))
RETURN
END FUNCTION iround



FUNCTION qf(alb, anc, n, irr, sigma, cc, lim1, acc, ith, trace, ifault)  &
   RESULT(fn_val)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2004-01-06  Time: 08:37:41

!     ALGORITHM AS 155  APPL. STATIST. (1980) VOL.29, NO.3

!     Distribution of a linear combination of non-central chi-squared
!     random variables.

IMPLICIT NONE
INTEGER, INTENT(IN)   :: irr
REAL, INTENT(IN)      :: alb(irr)
REAL, INTENT(IN)      :: anc(irr)
INTEGER, INTENT(IN)   :: n(irr)
REAL, INTENT(IN)      :: sigma
REAL, INTENT(IN)      :: cc
INTEGER, INTENT(IN)   :: lim1
REAL, INTENT(IN)      :: acc
INTEGER, INTENT(OUT)  :: ith(irr)
REAL, INTENT(OUT)     :: trace(7)
INTEGER, INTENT(OUT)  :: ifault
REAL                  :: fn_val

INTEGER  :: j, nj, nt, ntm
REAL     :: acc1, almx, xlim, xnt, xntm
REAL     :: utx, tausq, sd, aintv, aintv1, x, up, un, d1, d2, alj, ancj
REAL (dp) :: aintl, ersm
REAL     :: pi, aln28, sigsq, almax, almin, amean, c
INTEGER  :: icount, ir, lim
LOGICAL  :: ndtsrt, fail
REAL, PARAMETER  :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, four = 4.0, &
    sixtn = 16.0, fourp5 = 4.5, pt07 = 0.07, pt2 = 0.2, quart = 0.25,  &
    ten = 10.0, pt33 = 0.33, pt67 = 0.67, pt75 = 0.75, onept5 = 1.5,  &
    three = 3.0, onept1 = 1.1

!     Setting constants in COMMON.   ALN28 = ln(2) / 8.

pi = 3.14159265358979
aln28 = 0.0866

c = cc
ir = irr
lim = lim1
trace(1:7) = zero
ifault = 0
icount = 0
aintl = zero
ersm = zero
fn_val = -one
acc1 = acc
ndtsrt = .true.
fail = .false.

!     Find mean, sd, max & min of ALB.
!     Check that parameter values are valid.

xlim = lim
sigsq = sigma**2
sd = sigsq
almax = zero
almin = zero
amean = zero
j = 1

20 IF (.NOT.(j <= ir)) GO TO 60
nj = n(j)
alj = alb(j)
ancj = anc(j)
IF (.NOT.(nj < 0 .OR. ancj < zero)) GO TO 30
ifault = 3
GO TO 260

30 sd = sd + alj**2*(2*nj + four*ancj)
amean = amean + alj*(nj + ancj)
IF (.NOT.(almax < alj)) GO TO 40
almax = alj
GO TO 50

40 IF (.NOT.(almin > alj)) GO TO 50
almin = alj
50 j = j + 1
GO TO 20

60 IF (.NOT.(sd == zero)) GO TO 80
IF (.NOT.(c > zero)) GO TO 70
fn_val = one
GO TO 260

70 fn_val = zero
GO TO 260
80 IF (.NOT.(almin == zero .AND. almax == zero .AND. sigma == zero)) GO TO 90
ifault = 3
GO TO 260

90 sd = SQRT(sd)
IF (.NOT.(almax < -almin)) GO TO 100
almx = -almin
GO TO 110
100 almx = almax

!     Starting values for FINDU * CTFF.

110 utx = sixtn/sd
up = fourp5/sd
un = -up

!     Truncation point with no convergence factor.

CALL findu (n, alb, anc, utx, half*acc1)

!     Does convergence factor help ?

IF (.NOT.(c /= zero .AND. almx > pt07*sd)) GO TO 130
tausq = quart*acc1/cfe(n, alb, anc, ith, c)
IF (.NOT.(fail)) GO TO 120
fail = .false.
GO TO 130

120 IF (.NOT.(truncn(n, alb, anc, utx, tausq) < pt2*acc1)) GO TO 130
sigsq = sigsq + tausq
CALL findu (n, alb, anc, utx, quart*acc1)
trace(6) = SQRT(tausq)
130 trace(5) = utx
acc1 = half*acc1

!     Find 'range' of distribution, quit if outside of this.

140 d1 = ctff(n, alb, anc, acc1, up) - c
IF (.NOT.(d1 < zero)) GO TO 150
fn_val = one
GO TO 260
150 d2 = c - ctff(n, alb, anc, acc1, un)
IF (.NOT.(d2 < zero)) GO TO 160
fn_val = zero
GO TO 260

!     Find integration interval.

160 IF (.NOT.(d1 > d2)) GO TO 170
aintv = d1
GO TO 180
170 aintv = d2
180 aintv = two*pi/aintv

!     Calculate number of terms required for main & auxiliary
!     integrations.

xnt = utx/aintv
xntm = three/SQRT(acc1)
IF (.NOT.(xnt > xntm*onept5)) GO TO 220
IF (.NOT.(xntm > xlim)) GO TO 190
ifault = 1
GO TO 260

!     Parameters for auxiliary integration.

190 ntm = iround(xntm)
aintv1 = utx/xntm
x = two*pi/aintv1
IF (.NOT.(x <= ABS(c))) GO TO 200
GO TO 220

!     Calculate convergence factor.

200 tausq = cfe(n, alb, anc, ith, c - x) + cfe(n, alb, anc, ith, c + x)
tausq = pt33*acc1/(onept1*tausq)
IF (.NOT.(fail)) GO TO 210
GO TO 220
210 acc1 = pt67*acc1

!     Auxiliary integration.

CALL integr (n, alb, anc, ntm, aintv1, tausq, .false.)
xlim = xlim - xntm
sigsq = sigsq + tausq
trace(3) = trace(3) + 1
trace(2) = trace(2) + ntm + 1

!     Find truncation point with new convergence factor.

CALL findu (n, alb, anc, utx, quart*acc1)
acc1 = pt75*acc1
GO TO 140

!     Main integration.

220 trace(4) = aintv
IF (.NOT.(xnt > xlim)) GO TO 230
ifault = 1
GO TO 260

230 nt = iround(xnt)
CALL integr (n, alb, anc, nt, aintv, zero, .true.)
trace(3) = trace(3) + 1
trace(2) = trace(2) + nt + 1
fn_val = half - aintl
trace(1) = ersm
up = ersm

!     Test whether round-off error could be significant.
!     Allow for radix 8 or 16 machines.

x = up + acc/ten
j = 1
240 IF (.NOT.(j <= 8)) GO TO 260
IF (.NOT.(j*x == j*up)) GO TO 250
ifault = 2
250 j = j*2
GO TO 240

260 trace(7) = icount
RETURN
END FUNCTION qf



SUBROUTINE countr()

!     Count number of calls to ERRBD, TRUNCN & CFE.

IMPLICIT NONE
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim

icount = icount + 1
IF (.NOT.(icount > lim)) GO TO 20
WRITE (6, 10)
10 FORMAT (' qf: cannot locate integration parameters'/)
STOP
20 RETURN
END SUBROUTINE countr



FUNCTION alog1(x, first) RESULT(fn_val)

!     If FIRST then return ln(1 + x) else ln(1 + x) - x.

IMPLICIT NONE
REAL, INTENT(IN)     :: x
LOGICAL, INTENT(IN)  :: first
REAL                 :: fn_val

REAL    :: s, s1, term, y, ak
REAL, PARAMETER  :: pt1 = 0.1, one = 1.0, two = 2.0, three = 3.0

IF (.NOT.(ABS(x) > pt1)) GO TO 20
IF (.NOT.(first)) GO TO 10
fn_val = LOG(one + x)
GO TO 70

10 fn_val = LOG(one + x) - x
GO TO 70

20 y = x/(two + x)
term = two*y**3
ak = three
IF (.NOT.(first)) GO TO 30
s = two
GO TO 40

30 s = -x
40 s = s*y
y = y**2
s1 = s + term/ak

50 IF (.NOT.(s1 /= s)) GO TO 60
ak = ak + two
term = term*y
s = s1
s1 = s + term/ak
GO TO 50

60 fn_val = s
70 RETURN
END FUNCTION alog1



FUNCTION exp1(x) RESULT(fn_val)

IMPLICIT NONE
REAL, INTENT(IN)  :: x
REAL              :: fn_val

REAL, PARAMETER  :: zero = 0.0, neg50 = -50.0

IF (.NOT.(x < neg50)) GO TO 10
fn_val = zero
GO TO 20
10 fn_val = EXP(x)
20 RETURN
END FUNCTION exp1



SUBROUTINE order (alb, ith)

!     Find order of absolute values of ALB.

IMPLICIT NONE
REAL, INTENT(IN)         :: alb(*)
INTEGER, INTENT(IN OUT)  :: ith(*)

INTEGER  :: j, k, k1, ithk
REAL     :: alj
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim

j = 1
10 IF (.NOT.(j <= ir)) GO TO 60
alj = ABS(alb(j))
k = j - 1
20 IF (.NOT.(k > 0)) GO TO 40
ithk = ith(k)
k1 = k + 1
IF (.NOT.(alj > ABS(alb(ithk)))) GO TO 50
ith(k1) = ithk
k = k - 1
GO TO 20

40 k = 0
k1 = 1
50 ith(k1) = j
j = j + 1
GO TO 10

60 ndtsrt = .false.
RETURN
END SUBROUTINE order



FUNCTION errbd(n, alb, anc, uu, cx) RESULT(fn_val)

!     Find bound on tail probability using mgf.
!     Cut-off point returned to CX.

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n(*)
REAL, INTENT(IN)     :: alb(*)
REAL, INTENT(IN)     :: anc(*)
REAL, INTENT(IN)     :: uu
REAL, INTENT(OUT)    :: cx
REAL                 :: fn_val

REAL     :: u
REAL     :: sum1, alj, ancj, x, y, const
INTEGER  :: j, nj
! COMMON /qfcom/ aintl, ersm, pi, aln28, sigsq, almax, almin, amean, c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: half = 0.5, one = 1.0, two = 2.0

CALL countr()
u = uu
const = u*sigsq
sum1 = u*const
u = two*u
j = ir
10 IF (.NOT.(j > 0)) GO TO 20
nj = n(j)
alj = alb(j)
ancj = anc(j)
x = u*alj
y = one - x
const = const + alj*(ancj/y + nj)/y
sum1 = sum1 + ancj*(x/y)**2
sum1 = sum1 + nj*(x**2/y + alog1(-x,.false.))
j = j - 1
GO TO 10
20 fn_val = exp1(-half*sum1)
cx = const
RETURN
END FUNCTION errbd



FUNCTION ctff(n, alb, anc, accx, upn) RESULT(fn_val)

!     Find CTFF so that P(QF > CTFF) < ACCX if UPN > 0;
!     P(QF < CTFF) < ACCX otherwise.

IMPLICIT NONE
INTEGER, INTENT(IN)   :: n(*)
REAL, INTENT(IN)      :: alb(*)
REAL, INTENT(IN)      :: anc(*)
REAL, INTENT(IN)      :: accx
REAL, INTENT(IN OUT)  :: upn
REAL                  :: fn_val

REAL     :: u1, u2, u, rb, const, c1, c2
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: zero = 0.0, one = 1.0, two = 2.0

u2 = upn
u1 = zero
c1 = amean
IF (.NOT.(u2 > zero)) GO TO 10
rb = almax
GO TO 20
10 rb = almin
20 rb = two*rb
u = u2/(one + u2*rb)
30 IF (.NOT.(errbd(n, alb, anc, u, c2) > accx)) GO TO 40
u1 = u2
c1 = c2
u2 = two*u2
u = u2/(one + u2*rb)
GO TO 30

40 u = (c1 - amean)/(c2 - amean)
50 IF (.NOT.(u < 0.9)) GO TO 80
u = (u1 + u2)/two
IF (.NOT.(errbd(n, alb, anc, u/(one + u*rb), const) > accx)) GO TO 60
u1 = u
c1 = const
GO TO 70

60 u2 = u
c2 = const
70 u = (c1 - amean)/(c2 - amean)
GO TO 50

80 fn_val = c2
upn = u2
RETURN
END FUNCTION ctff



FUNCTION truncn(n, alb, anc, uu, tausq) RESULT(fn_val)

!     Bound integration error due to truncation at U.

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n(*)
REAL, INTENT(IN)     :: alb(*)
REAL, INTENT(IN)     :: anc(*)
REAL, INTENT(IN)     :: uu
REAL, INTENT(IN)     :: tausq
REAL                 :: fn_val

REAL     :: u
REAL     :: sum1, sum2, prod1, prod2, prod3, alj, ancj, x, y, err1, err2
INTEGER  :: j, nj, ns
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: zero = 0.0, quart = 0.25, half = 0.5, one = 1.0, two = 2.0

CALL countr()
u = uu
sum1 = zero
prod2 = zero
prod3 = zero
ns = 0
sum2 = (sigsq + tausq)*u**2
prod1 = two*sum2
u = two*u
j = 1

10 IF (.NOT.(j <= ir)) GO TO 40
alj = alb(j)
ancj = anc(j)
nj = n(j)
x = (u*alj)**2
sum1 = sum1 + ancj*x/(one + x)
IF (.NOT.(x > one)) GO TO 20
prod2 = prod2 + nj*LOG(x)
prod3 = prod3 + nj*alog1(x,.true.)
ns = ns + nj
GO TO 30

20 prod1 = prod1 + nj*alog1(x,.true.)
30 j = j + 1
GO TO 10

40 sum1 = half*sum1
prod2 = prod1 + prod2
prod3 = prod1 + prod3
x = exp1(-sum1 - quart*prod2)/pi
y = exp1(-sum1 - quart*prod3)/pi
IF (.NOT.(ns == 0)) GO TO 50
err1 = one
GO TO 60

50 err1 = x*two/ns
60 IF (.NOT.(prod3 > one)) GO TO 70
err2 = 2.5*y
GO TO 80

70 err2 = one
80 IF (.NOT.(err2 < err1)) GO TO 90
err1 = err2
90 x = half*sum2
IF (.NOT.(x <= y)) GO TO 100
err2 = one
GO TO 110

100 err2 = y/x
110 IF (.NOT.(err1 < err2)) GO TO 120
fn_val = err1
GO TO 130

120 fn_val = err2
130 RETURN
END FUNCTION truncn



SUBROUTINE findu (n, alb, anc, utx, accx)

!     Find U such that TRUNCN(U) < ACCX & TRUNCN(U / 1.2) > ACCX.

IMPLICIT NONE
INTEGER, INTENT(IN)   :: n(*)
REAL, INTENT(IN)      :: alb(*)
REAL, INTENT(IN)      :: anc(*)
REAL, INTENT(IN OUT)  :: utx
REAL, INTENT(IN)      :: accx

REAL     :: u, ut
INTEGER  :: i
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: divis(4) = (/ 2.0, 1.4, 1.2, 1.1 /)
REAL, PARAMETER  :: four = 4.0, zero = 0.0

ut = utx
u = ut/four
IF (.NOT.(truncn(n, alb, anc, u, zero) > accx)) GO TO 20
u = ut
10 IF (.NOT.(truncn(n, alb, anc, u, zero) > accx)) GO TO 40
ut = ut*four
u = ut
GO TO 10

20 ut = u
u = u/four
30 IF (.NOT.(truncn(n, alb, anc, u, zero) <= accx)) GO TO 40
ut = u
u = u/four
GO TO 30

40 DO  i = 1,4
  u = ut/divis(i)
  IF (.NOT.(truncn(n, alb, anc, u, zero) <= accx)) CYCLE
  ut = u
END DO
utx = ut
RETURN
END SUBROUTINE findu



SUBROUTINE integr (n, alb, anc, nterm, aintrv, tausq, main)

!     Carry out integration with NTERM terms, at interval AINTRV.
!     If not MAIN then multiply integrand by  1 - exp(-0.5 * TAUSQ * U**2).

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n(*)
REAL, INTENT(IN)     :: alb(*)
REAL, INTENT(IN)     :: anc(*)
INTEGER, INTENT(IN)  :: nterm
REAL, INTENT(IN)     :: aintrv
REAL, INTENT(IN)     :: tausq
LOGICAL, INTENT(IN)  :: main

REAL      :: ainpi, u, sum1, sum2, sum3, x, y, z
INTEGER   :: k, j, nj
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: quart = 0.25, half = 0.5, one = 1.0, two = 2.0

ainpi = aintrv/pi
k = nterm
10 IF (.NOT.(k >= 0)) GO TO 50
u = (k + half)*aintrv
sum1 = -two*u*c
sum2 = ABS(sum1)
sum3 = -half*sigsq*u**2
j = ir

20 IF (.NOT.(j > 0)) GO TO 30
nj = n(j)
x = two*alb(j)*u
y = x**2
sum3 = sum3 - quart*nj*alog1(y,.true.)
y = anc(j)*x/(one + y)
z = nj*ATAN(x) + y
sum1 = sum1 + z
sum2 = sum2 + ABS(z)
sum3 = sum3 - half*x*y
j = j - 1
GO TO 20

30 x = ainpi*exp1(sum3)/u
IF (.NOT.(.NOT.main)) GO TO 40
x = x*(one - exp1(-half*tausq*u**2))
40 sum1 = SIN(half*sum1)*x
sum2 = half*sum2*x
aintl = aintl + sum1
ersm = ersm + sum2
k = k - 1
GO TO 10

50 RETURN
END SUBROUTINE integr



FUNCTION cfe(n, alb, anc, ith, x) RESULT(fn_val)

!     Coefficient of TAUSQ in error when convergence factor of
!     exp(-0.5 * TAUSQ * U**2) is used when df is evaluated at X.

IMPLICIT NONE
INTEGER, INTENT(IN)   :: n(*)
REAL, INTENT(IN)      :: alb(*)
REAL, INTENT(IN)      :: anc(*)
INTEGER, INTENT(OUT)  :: ith(*)
REAL, INTENT(IN)      :: x
REAL                  :: fn_val

REAL      :: axl, axl1, axl2, sxl, sum1, alj
INTEGER   :: j, k, it, itk
! COMMON /qfcom/ aintl,ersm,pi,aln28,sigsq,almax,almin,amean,c,  &
!    icount,ir,ndtsrt,fail,lim
REAL, PARAMETER  :: zero = 0.0, one = 1.0, four = 4.0, hundrd = 100.0

CALL countr()
IF (.NOT.(ndtsrt)) GO TO 10
CALL order (alb, ith)
10 axl = ABS(x)
sxl = SIGN(one,x)
sum1 = zero
j = ir

20 IF (.NOT.(j > 0)) GO TO 70
it = ith(j)
IF (.NOT.(alb(it)*sxl > zero)) GO TO 60
alj = ABS(alb(it))
axl1 = axl - alj*(n(it) + anc(it))
axl2 = alj/aln28
IF (.NOT.(axl1 > axl2)) GO TO 30
axl = axl1
GO TO 60

30 IF (.NOT.(axl > axl2)) GO TO 40
axl = axl2
40 sum1 = (axl - axl1)/alj
k = j - 1
50 IF (.NOT.(k > 0)) GO TO 70
itk = ith(k)
sum1 = sum1 + (n(itk) + anc(itk))
k = k - 1
GO TO 50

60 j = j - 1
GO TO 20

70 IF (.NOT.(sum1 > hundrd)) GO TO 80
fn_val = one
fail = .true.
GO TO 90
80 fn_val = 2**(sum1/four)/(pi*axl**2)
90 RETURN
END FUNCTION cfe

END MODULE Nonc_chisq