SUBROUTINE FUNPRE (FUNC,XFCUM,X2LOW,X2HIGH)
C         F. JAMES,    MAY, 1976
C         MODIFIED OCT, 1980 TO ADD PRINTOUT OF INTEGRAL
C         MODIFIED DEC., 1980 TO DELETE LEADING AND TRAILING
C            RANGES OF X WHERE FUNCTION IS ZERO.
C         MODIFIED JUNE,1982 TO FIX POSSIBLE INFINITE LOOP.
C
C         PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN
C         BY FINDING THE PERCENTILES
C         (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION)
      EXTERNAL FUNC
      COMMON/FUNINT/TFTOT
      DIMENSION XFCUM(100)
      DATA NBINS / 99/
      DATA NZ / 10/
      DATA MAXZ / 20/
      DATA NITMAX / 6 /
      DATA IFUNC/0/
      IFUNC = IFUNC + 1
C         FIND MACHINE ACCURACY
      COMP1 = 1.0
      DO 200 I= 1, 100
      COMP1 = COMP1*0.5
      COMP2 = 1.0 - COMP1
      IF(COMP2 .EQ. 1.0) GOTO 210
  200 CONTINUE
      COMP1 = 1.0E-10
  210 PRECIS = COMP1
C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
      CALL FUNZER(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
      XRANGE = XHIGH-XLOW
      IF(XRANGE .LE. 0.) GOTO 900
      RTEPS = MAX(0.0001,PRECIS*10.)
      TFTOT = GAUSS(FUNC,XLOW,XHIGH,RTEPS)
C         PRINT OUT VALUE OF NORMALIZATION INTEGRAL
      WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
      RTEPS = 0.001
      IF(TFTOT .LE. 0.) GOTO 900
      TPCTIL = TFTOT/NBINS
      TZ = TPCTIL/NZ
      TZMAX = TZ * 2.
       XFCUM(1) = XLOW
      XFCUM(NBINS+1) = XHIGH
      X = XLOW
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
      NBINM1 = NBINS - 1
C         LOOP OVER BINS (HUNDREDTH PERCENTILES)
      DO 600 IBIN = 1, NBINM1
      TCUM = 0.
      X1 = X
      F1 = F
      DXMAX = (XHIGH -X) / NZ
      FMIN = TZ/DXMAX
      FMINZ = FMIN
C         LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE
      DO 500 IZ= 1, MAXZ
      XINCR = TZ/MAX(F1,FMIN,FMINZ)
  350 X = X1 + XINCR
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
      TINCR = (X-X1) * 0.5 * (F+F1)
      IF(TINCR .LT. TZMAX) GOTO 370
      XINCR = XINCR * 0.5
      GOTO 350
  370 CONTINUE
      TCUM = TCUM + TINCR
      IF(TCUM .GE. TPCTIL*0.99) GOTO 520
      FMINZ = TZ*F/ (TPCTIL-TCUM)
      F1 = F
      X1 = X
  500 CONTINUE
      WRITE(6,2000)
 2000 FORMAT('0      FAILURE TO FIND TRAPEZOID   HELP')
C         END OF TRAPEZOID LOOP
C         ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION
  520 CONTINUE
      X1 = XFCUM(IBIN)
      XBEST = X
      DTBEST = TPCTIL
      TPART = TPCTIL
C         ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS
      DO 550 IHOME= 1, NITMAX
      XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
  535 X = XBEST + XINCR
      X2 = X
      TPART2 = GAUSS(FUNC,X1,X2,RTEPS)
      DTPAR2 = TPART2-TPCTIL
      DTABS = ABS(DTPAR2)
      IF(ABS(XINCR) .LT. PRECIS) GOTO 545
      IF(DTABS .LT. DTBEST) GOTO 545
      XINCR = XINCR * 0.5
      GOTO 535
  545 DTBEST = DTABS
      XBEST = X
      IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
      TPART = TPART2
      F = FUNC(X)
  550 CONTINUE
      IHOME = NITMAX
C
  580 CONTINUE
      XFCUM(IBIN+1) = X
      F = FUNC(X)
      IF(F .LT. 0.) GOTO 900
  600 CONTINUE
C         END OF LOOP OVER BINS
      X1 = XFCUM(NBINS)
      X2 = XHIGH
      TPART = GAUSS(FUNC,X1,X2,RTEPS)
      ABERR = ABS(TPART-TPCTIL)/TFTOT
      WRITE(6,1001) IFUNC,ABERR
      IF(ABERR .GT. RTEPS)  WRITE(6,1002)
      RETURN
  900 WRITE(6,1000) X,F,XLOW,XHIGH
      RETURN
 1000 FORMAT('0FUNPRE FINDS NEGATIVE FUNCTION VALUE OR RANGE OF X'/
     + ,' X=',E15.6,', F=',E15.6,20X,'XLOW=',E15.6,' XHIGH=',E15.6/)
 1001 FORMAT(' SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION NUMBER',I4,
     + ' FOR FUNRAN'/' MAXIMUM RELATIVE ERROR IN CUMULATIVE',
     + 'DISTRIBUTION WILL BE',E15.5)
 1002 FORMAT('+',80X,'WARNING,THIS MAY BE TOO BIG'//)
 1003 FORMAT('0SUBROUTINE FUNPRE FINDS THE INTEGRAL OF USER FUNCTION',
     +  I2,' FROM ',E12.5,' TO ',E12.5,' IS ',E14.6)
      END