*     G05YBF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO=0.0D0)
*     .. Local Scalars ..
      DOUBLE PRECISION SUM, VAL1, VAL2
      INTEGER          I, IDIM, IFAIL, NTIMES, SKIP
      LOGICAL          LNORM
*     .. Local Arrays ..
      DOUBLE PRECISION MEAN(20), QUASI(20), STD(20)
      INTEGER          IREF(2000)
*     .. External Functions ..
      DOUBLE PRECISION FUN
      EXTERNAL         FUN
*     .. External Subroutines ..
      EXTERNAL         G05YBF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'G05YBF Example Program Results'
      NTIMES = 10000
      IDIM = 20
      LNORM = .FALSE.
      DO 20 I = 1, IDIM
         MEAN(I) = 2.0D0
         STD(I) = 1.0D0
   20 CONTINUE
*
      CALL G05YBF(.TRUE.,'F',LNORM,MEAN,STD,SKIP,IDIM,QUASI,IREF,IFAIL)
*
      SUM = ZERO
      DO 40 I = 1, NTIMES
         IFAIL = 0
*
         CALL G05YBF(.FALSE.,'F',LNORM,MEAN,STD,SKIP,IDIM,QUASI,IREF,
     +               IFAIL)
*
         SUM = SUM + FUN(IDIM,MEAN,STD,QUASI)
   40 CONTINUE
      VAL1 = SUM/DBLE(NTIMES)
      WRITE (NOUT,*)
      WRITE (NOUT,99999) 'Calculate value of the integral = ', VAL1
      VAL2 = DBLE(IDIM)*3.0D0
      WRITE (NOUT,*)
      WRITE (NOUT,99999) 'Exact value of the integral     = ', VAL2
      STOP
*
99999 FORMAT (1X,A,F8.3)
      END
*
      DOUBLE PRECISION FUNCTION FUN(IDIM,MEAN,STD,X)
*     .. Parameters ..
      DOUBLE PRECISION              ZERO
      PARAMETER                     (ZERO=0.0D0)
*     .. Scalar Arguments ..
      INTEGER                       IDIM
*     .. Array Arguments ..
      DOUBLE PRECISION              MEAN(IDIM), STD(IDIM), X(IDIM)
*     .. Local Scalars ..
      DOUBLE PRECISION              TMP1, TMP2
      INTEGER                       J
*     .. Executable Statements ..
      TMP1 = ZERO
      DO 20 J = 1, IDIM
         TMP2 = (X(J)-MEAN(J))/STD(J)
         TMP1 = TMP1 + TMP2*TMP2*TMP2*TMP2
   20 CONTINUE
      FUN = TMP1
      RETURN

      END