REAL*8 FUNCTION YUM(F,C1,C2,XSTAR,D2LGF)
C  This subroutine draws from a distribution with a unimodal density on a
C  finite support.  Documentation is in notes on the ARMA(p,q) paper.
C  Inputs:
C  F      Function F(X) that has input ordinate X and returns density kernel F.
C  C1     Lower end of finite support
C  C2     Upper end of finite support
C  XSTAR  Mode of density kernel
C  D2LGF  Second derivative of log density kernel at mode
C  Output:
C  YUM    Drawing from distribution
      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL F
      DATA TSTAR,E1,E2,IR/10.0D0,0.0D0,1.0D-8,2/
C  Step 1
  100 CALL DQDAG(F,C1,XSTAR,E1,E2,IR,R1,E3)
      CALL DQDAG(F,XSTAR,C2,E1,E2,IR,R2,E3)
      R=R1+R2
      A=1.0D0/DSQRT(-D2LGF)
      B1P=XSTAR
      B2M=XSTAR
C  Step 2
  200 SSTAR=R*DRNUNF()
      S=0.0D0
C  Step 3
  300 IF(B1P.GE.C2)GO TO 400
  310 Q=A
      F1P=F(B1P)
      QMAX=C2-B1P
      IF(Q.GT.QMAX)Q=QMAX
  320 B2P=B1P+Q
      IF((TSTAR*F(B2P)).GT.F1P)GO TO 330
      Q=.5D0*Q
      GO TO 320
  330 CALL DQDAG(F,B1P,B2P,E1,E2,IR,SINC,E3)
      S=S+SINC
      IF(S.GE.SSTAR)GO TO 340
      B1P=B2P
      GO TO 400
  340 T=F1P
      B1=B1P
      B2=B2P
      GO TO 500
C  Step 4
  400 IF(B2M.GT.C1)GO TO 410
      IF(B1P.GE.C2)GO TO 610
      GO TO 310
  410 Q=A
      F2M=F(B2M)
      QMAX=B2M-C1
      IF(Q.GT.QMAX)Q=QMAX
  420 B1M=B2M-Q
      IF((TSTAR*F(B1M)).GT.F2M)GO TO 430
      Q=.5D0*Q
      GO TO 420
  430 CALL DQDAG(F,B1M,B2M,E1,E2,IR,SINC,E3)
      S=S+SINC
      IF(S.GT.SSTAR)GO TO 440
      B2M=B1M
      GO TO 300
  440 T=F2M
      B1=B1M
      B2=B2M
C  Step 5
  500 X=B1+(B2-B1)*DRNUNF()
      IF(T*DRNUNF().GT.F(X))GO TO 500
      YUM=X
      RETURN
C  Errors
  610 WRITE(6,620)
  620 FORMAT(/,'  Adding-up violations in YUM.')
      STOP
      END