C      ALGORITHM 451, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 16, NO. 8, August, 1973, PP.483--485.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Thu Dec 15 13:28:14 2005
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver
RESULTS= Res

Objs1= driver.o src.o
driver: $(Objs1)
	$(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS)
Res: driver 
	./driver >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c***********************************************************************
c
cc TOMS451_PRB tests CHISQD.
c
      implicit none

      integer n_num
      integer p_num

      parameter ( n_num = 10 )
      parameter ( p_num = 5 )

      real chisqd
      integer i
      integer j
      integer n
      integer n_vec(n_num)
      real p
      real p_vec(p_num)
      real row(p_num)
      real table(n_num,p_num)

      save n_vec
      save p_vec
      save table

      data n_vec / 1, 2, 3, 4, 5, 10, 15, 20, 50, 100 /
      data p_vec / 
     &  0.9995E+00, 0.9950E+00, 0.500E+00, 0.0010E+00, 0.0001E+00 /
c
c  The data for TABLE is listed one column at a time.  That is,
c  a consecutive pair of rows of data represent a COLUMN of the array.
c
      data table / 
     &   0.000000,  0.001000,  0.015312,  0.063955,   0.158168,
     &   1.264941,  3.107881,  5.398208, 23.460876,  59.895508,
     &   0.000039,  0.010025,  0.071641,  0.206904,   0.411690,
     &   2.155869,  4.601008,  7.433892, 27.990784,  67.327621,
     &   0.454933,  1.386293,  2.365390,  3.356400,   4.351295,
     &   9.341794, 14.338853, 19.337418, 49.334930,  99.334122,
     &  10.827576, 13.815512, 16.268982, 18.467987,  20.515503,
     &  39.589081, 37.697662, 45.314896, 86.660767, 149.449051,
     &  15.135827, 18.420670, 21.106873, 23.510040,  25.744583,
     &  35.565170, 44.267853, 52.387360, 95.969482, 161.319733 /

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS451_PRB'
      write ( *, '(a)' ) '  Test TOMS algorithm 451, which computes'
      write ( *, '(a)' ) '  values of the Chi-square quantile,'
      write ( *, '(a)' ) '    CHISQD ( P, N )'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  We print out a table using both CHISQD'
      write ( *, '(a)' ) '  and the values printed in the writeup.'
      write ( *, '(a)' ) ' '
      write ( *, '(5x,5(2x,f12.6))' ) ( p_vec(j), j = 1, p_num )
      write ( *, '(a)' ) ' '

      do i = 1, n_num

        n = n_vec(i)

        do j = 1, p_num

          p = p_vec(j)
          row(j) = chisqd ( p, n )

        end do

        write ( *, '(2x,i3,5(2x,f12.6))' ) 
     &    n_vec(i), ( row(j), j = 1, p_num )
        write ( *, '(2x,i3,5(2x,f12.6))' ) 
     &    n_vec(i), ( table(i,j), j = 1, p_num )

      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS451_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'

      stop
      end
      function gaussd ( p )

!***********************************************************************
!
!! GAUSSD inverts the standard normal CDF.
!
!  Discussion:
!
!    The result is accurate to about 1 part in 10**16.
!
!  Modified:
!
!    27 December 2004
!
!  Author:
!
!    Michael Wichura
!
!  Reference:
!
!    Michael Wichura,
!    The Percentage Points of the Normal Distribution,
!    Algorithm AS 241,
!    Applied Statistics,
!    Volume 37, Number 3, pages 477-484, 1988.
!
!  Parameters:
!
!    Input, real P, the value of the cumulative probability
!    densitity function.  0 < P < 1.  If P is outside this range, an
!    "infinite" value will be returned.
!
!    Output, real GAUSSD, the normal deviate value
!    with the property that the probability of a standard normal deviate being
!    less than or equal to the value is P.
!
      implicit none

      real a(8)
      real b(8)
      real c(8)
      real const1
      real const2
      real d(8)
      real e(8)
      real f(8)
      real gaussd
      real p
      real q
      real r
      real rpoly_value
      real split1
      real split2

      save a
      save b
      save c
      save const1
      save const2
      save d
      save e
      save f
      save split1
      save split2

      data a /
     &  3.3871328727963666080E+00, 
     &  1.3314166789178437745E+02, 
     &  1.9715909503065514427E+03, 
     &  1.3731693765509461125E+04, 
     &  4.5921953931549871457E+04, 
     &  6.7265770927008700853E+04, 
     &  3.3430575583588128105E+04, 
     &  2.5090809287301226727E+03 /
      data b /
     &  1.0E+00, 
     &  4.2313330701600911252E+01, 
     &  6.8718700749205790830E+02, 
     &  5.3941960214247511077E+03, 
     &  2.1213794301586595867E+04, 
     &  3.9307895800092710610E+04, 
     &  2.8729085735721942674E+04, 
     &  5.2264952788528545610E+03 /
      data c /
     &  1.42343711074968357734E+00, 
     &  4.63033784615654529590E+00, 
     &  5.76949722146069140550E+00, 
     &  3.64784832476320460504E+00, 
     &  1.27045825245236838258E+00, 
     &  2.41780725177450611770E-01, 
     &  2.27238449892691845833E-02, 
     &  7.74545014278341407640E-04 /
      data const1 / 0.180625E+00 /
      data const2 / 1.6E+00 /
      data d /
     &  1.0E+00, 
     &  2.05319162663775882187E+00, 
     &  1.67638483018380384940E+00, 
     &  6.89767334985100004550E-01, 
     &  1.48103976427480074590E-01, 
     &  1.51986665636164571966E-02, 
     &  5.47593808499534494600E-04, 
     &  1.05075007164441684324E-09 /
      data e /
     &  6.65790464350110377720E+00, 
     &  5.46378491116411436990E+00, 
     &  1.78482653991729133580E+00, 
     &  2.96560571828504891230E-01, 
     &  2.65321895265761230930E-02, 
     &  1.24266094738807843860E-03, 
     &  2.71155556874348757815E-05, 
     &  2.01033439929228813265E-07 /
      data f /
     &  1.0E+00, 
     &  5.99832206555887937690E-01, 
     &  1.36929880922735805310E-01, 
     &  1.48753612908506148525E-02, 
     &  7.86869131145613259100E-04, 
     &  1.84631831751005468180E-05, 
     &  1.42151175831644588870E-07, 
     &  2.04426310338993978564E-15 /
      data split1 / 0.425E+00 /
      data split2 / 5.0E+00 /

      if ( p <= 0.0E+00 ) then
        gaussd = -1.0E+30
        return
      end if

      if ( 1.0D+00 <= p ) then
        gaussd = 1.0E+30
        return
      end if

      q = p - 0.5E+00

      if ( abs ( q ) <= split1 ) then

        r = const1 - q * q
        gaussd = q * rpoly_value ( 8, a, r ) / rpoly_value ( 8, b, r )

      else

        if ( q < 0.0E+00 ) then
          r = p
        else
          r = 1.0E+00 - p
        end if

        if ( r <= 0.0E+00 ) then
          gaussd = -1.0E+00
          stop
        end if

        r = sqrt ( -log ( r ) )

        if ( r <= split2 ) then

          r = r - const2
          gaussd = rpoly_value ( 8, c, r ) / rpoly_value ( 8, d, r )

        else

          r = r - split2
          gaussd = rpoly_value ( 8, e, r ) / rpoly_value ( 8, f, r )

        end if

        if ( q < 0.0E+00 ) then
          gaussd = -gaussd
        end if

      end if

      return
      end
      function rpoly_value ( n, a, x )

c***********************************************************************
c
cc RPOLY_VALUE evaluates a double precision polynomial.
c
c  Discussion:
c
c    For sanity's sake, the value of N indicates the NUMBER of
c    coefficients, or more precisely, the ORDER of the polynomial,
c    rather than the DEGREE of the polynomial.  The two quantities
c    differ by 1, but cause a great deal of confusion.
c
c    Given N and A, the form of the polynomial is:
c
c      p(x) = a(1) + a(2) * x + ... + a(n-1) * x^(n-2) + a(n) * x^(n-1)
c
c  Modified:
c
c    13 August 2004
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer N, the order of the polynomial.
c
c    Input, real A(N), the coefficients of the polynomial.
c    A(1) is the constant term.
c
c    Input, real X, the point at which the polynomial is
c    to be evaluated.
c
c    Output, real RPOLY_VALUE, the value of the polynomial at X.
c
      implicit none

      integer n

      real a(n)
      integer i
      real rpoly_value
      real x

      rpoly_value = 0.0E+00
      do i = n, 1, -1
        rpoly_value = rpoly_value * x + a(i)
      end do

      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS451_PRB
  Test TOMS algorithm 451, which computes
  values of the Chi-square quantile,
    CHISQD ( P, N )
 
  We print out a table using both CHISQD
  and the values printed in the writeup.
 
           0.999500      0.995000      0.500000      0.001000      0.000100
 
    1      0.000000      0.000039      0.454936     10.827565     15.136707
    1      0.000000      0.000039      0.454933     10.827576     15.135827
    2      0.001000      0.010025      1.386294     13.815511     18.420681
    2      0.001000      0.010025      1.386293     13.815512     18.420670
    3      0.015313      0.071641      2.365390     16.269005     21.106068
    3      0.015312      0.071641      2.365390     16.268982     21.106873
    4      0.063959      0.206904      3.356400     18.468008     23.509207
    4      0.063955      0.206904      3.356400     18.467987     23.510040
    5      0.158176      0.411690      4.351294     20.515549     25.743729
    5      0.158168      0.411690      4.351295     20.515503     25.744583
   10      1.264977      2.155873      9.341793     29.589115     35.564262
   10      1.264941      2.155869      9.341794     39.589081     35.565170
   15      3.107936      4.601000     14.338851     37.697678     44.266891
   15      3.107881      4.601008     14.338853     37.697662     44.267853
   20      5.398283      7.433881     19.337423     45.314938     52.386230
   20      5.398208      7.433892     19.337418     45.314896     52.387360
   50     23.461065     27.990749     49.334942     86.660881     95.968155
   50     23.460876     27.990784     49.334930     86.660767     95.969482
  100     59.895790     67.327560     99.334129    149.449326    161.317841
  100     59.895508     67.327621     99.334122    149.449051    161.319733
 
TOMS451_PRB
  Normal end of execution.
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      REAL FUNCTION CHISQD(P, N)
      REAL C, A, P, F, F1, F2, T, GAUSSD
      INTEGER N
      DIMENSION C(21), A(19)
      DATA C(1)/1.565326E-3/, C(2)/1.060438E-3/,
     &  C(3)/-6.950356E-3/, C(4)/-1.323293E-2/,
     &  C(5)/2.277679E-2/, C(6)/-8.986007E-3/,
     &  C(7)/-1.513904E-2/, C(8)/2.530010E-3/,
     &  C(9)/-1.450117E-3/, C(10)/5.169654E-3/,
     &  C(11)/-1.153761E-2/, C(12)/1.128186E-2/,
     &  C(13)/2.607083E-2/, C(14)/-0.2237368/,
     &  C(15)/9.780499E-5/, C(16)/-8.426812E-4/,
     &  C(17)/3.125580E-3/, C(18)/-8.553069E-3/,
     &  C(19)/1.348028E-4/, C(20)/0.4713941/, C(21)/1.0000886/
      DATA A(1)/1.264616E-2/, A(2)/-1.425296E-2/,
     &  A(3)/1.400483E-2/, A(4)/-5.886090E-3/,
     &  A(5)/-1.091214E-2/, A(6)/-2.304527E-2/,
     &  A(7)/3.135411E-3/, A(8)/-2.728484E-4/,
     &  A(9)/-9.699681E-3/, A(10)/1.316872E-2/,
     &  A(11)/2.618914E-2/, A(12)/-0.2222222/,
     &  A(13)/5.406674E-5/, A(14)/3.483789E-5/,
     &  A(15)/-7.274761E-4/, A(16)/3.292181E-3/,
     &  A(17)/-8.729713E-3/, A(18)/0.4714045/, A(19)/1./
      IF (N-2) 10, 20, 30
10    CHISQD = GAUSSD(.5*P)
      CHISQD = CHISQD*CHISQD
      RETURN
20    CHISQD = -2. * ALOG(P)
      RETURN
30    F = N
      F1 = 1. / F
      T = GAUSSD(1.-P)
      F2 = SQRT(F1) * T
      IF ( N.GE.(2+INT(4.*ABS(T)))) GO TO 40
      CHISQD=(((((((C(1)*F2+C(2))*F2+C(3))*F2+C(4))*F2
     &  +C(5))*F2+C(6))*F2+C(7))*F1+((((((C(8)+C(9)*F2)*F2
     &  +C(10))*F2+C(11))*F2+C(12))*F2+C(13))*F2+C(14)))*F1 +
     &  (((((C(15)*F2+C(16))*F2+C(17))*F2+C(18))*F2
     &  +C(19))*F2+C(20))*F2+C(21)
      GO TO 50
40    CHISQD = (((A(1)+A(2)*F2)*F1+(((A(3)+A(4)*F2)*F2
     &  +A(5))*F2+A(6)))*F1+(((((A(7)+A(8)*F2)*F2+A(9))*F2
     &  +A(10))*F2+A(11))*F2+A(12)))*F1 + (((((A(13)*F2
     &  +A(14))*F2+A(15))*F2+A(16))*F2+A(17))*F2*F2
     &  +A(18))*F2+A(19)
50    CHISQD = CHISQD*CHISQD*CHISQD*F
      RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0