C $$                   ******  TEST1K  ******
C TEST PROGRAM FOR MP PACKAGE
C
C THIS PROGRAM COMPUTES THE CONSTANTS GIVEN IN APPENDIX A
C OF KNUTH, THE ART OF COMPUTER PROGRAMMING, VOL. 3.
C THE CONSTANTS ARE PRINTED IN THE SAME ORDER AS IN KNUTH.
C
C THE CONSTANTS ARE COMPUTED TO 1000 DECIMAL PLACES, BUT
C TO INCREASE THE ACCURACY IT IS ONLY NECESSARY TO CHANGE
C THE STATEMENT PLACES = 1000, AND POSSIBLY
C THE PARAMETERS OF THE CALL TO MPSET2 AND THE DIMENSIONS
C OF THE ARRAYS.
C
C THE MP ROUTINES NEEDED TO RUN THIS PROGRAM ARE LISTED IN THE MP
C USERS GUIDE.
C
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY
      INTEGER B, DUMMY(16), LUN, M, MXR, MXSPTR, SPTR, T
C DIMENSIONS OF R, X, ETC. CAN BE REDUCED IF WORDLENGTH .GT. 16 BITS.
      INTEGER R(3300)
C TEMPORARY MP VARIABLES REQUIRE SPACE T+2 AND T .LE. 248.
      INTEGER I, PLACES, X(250), Y(250), PHI(250), PI(250)
C
C SET OUTPUT UNIT = 6 AND WORKING PRECISION TO THE
C EQUIVALENT OF AT LEAST PLACES DECIMAL PLACES.  THE OTHER
C PARAMETERS ARE THE DIMENSIONS OF X AND THE INDICES OF THE
C FIRST AND LAST WORDS OF BLANK COMMON AVAILABLE TO MP.
      PLACES = 1000
      CALL MPSET2 (6, PLACES+3, 250, 1, 3300)
      WRITE (LUN, 5) B, T
    5 FORMAT (29H1TEST OF MP PACKAGE,   BASE =, I9,
     $        11H,  DIGITS =, I4 ///)
C
C COMPUTE SQRT(2), SQRT(3), SQRT(5) AND SQRT(10)
      DO 10 I = 2, 5
      CALL MPQPWR ((5*I)/4 + 4*(I/5), 1, 1, 2, X)
   10 CALL MP40D (PLACES, X)
C COMPUTE 2**(1/3) AND 3**(1/3)
      DO 20 I = 2, 3
      CALL MPQPWR (IABS(I), 1, 1, 3, X)
   20 CALL MP40D (PLACES, X)
C COMPUTE 2**(1/4)
      CALL MPQPWR (2, 1, 1, 4, X)
      CALL MP40D (PLACES, X)
C COMPUTE LN(2), LN(3) AND LN(10)
      DO 30 I = 2, 4
      CALL MPLNI (I + 6*(I/4), X)
   30 CALL MP40D (PLACES, X)
C COMPUTE 1/LN(2) AND 1/LN(10)
C COULD HAVE SAVED ABOVE RESULTS TO SPEED UP HERE
      DO 40 I = 1, 2
      CALL MPLNI (8*I - 6, X)
      CALL MPREC (X, X)
   40 CALL MP40D (PLACES, X)
C COMPUTE PI, PI/180, 1/PI, PI**2, SQRT(PI)
      CALL MPPI (PI)
      CALL MP40D (PLACES, PI)
      CALL MPDIVI (PI, 180, Y)
      CALL MP40D (PLACES, Y)
      CALL MPREC (PI, Y)
      CALL MP40D (PLACES, Y)
      CALL MPPWR (PI, 2, Y)
      CALL MP40D (PLACES, Y)
      CALL MPSQRT (PI, Y)
      CALL MP40D (PLACES, Y)
C COMPUTE GAMMA (1/3)
      CALL MPGAMQ (1, 3, X)
      CALL MP40D (PLACES, X)
C COMPUTE GAMMA (2/3) FROM GAMMA (1/3) (WE COULD
C ALSO CALL MPGAMQ (2, 3, X))
      CALL MPQPWR (3, 4, 1, 2, Y)
      CALL MPMUL (X, Y, X)
      CALL MPDIV (PI, X, X)
      CALL MP40D (PLACES, X)
C COMPUTE E, 1/E, AND E**2
      CALL MPCIM (1, X)
      CALL MPEXP (X, X)
      CALL MP40D (PLACES, X)
      CALL MPREC (X, Y)
      CALL MP40D (PLACES, Y)
      CALL MPMUL (X, X, Y)
      CALL MP40D (PLACES, Y)
C COMPUTE EULERS CONSTANT (GAMMA)
      CALL MPEUL (X)
      CALL MP40D (PLACES, X)
C COMPUTE LN(PI), PHI
      CALL MPLN (PI, Y)
      CALL MP40D (PLACES, Y)
      CALL MPQPWR (5, 4, 1, 2, Y)
      CALL MPADDQ (Y, 1, 2, PHI)
      CALL MP40D (PLACES, PHI)
C COMPUTE EXP(GAMMA)  (GAMMA IS IN X)
      CALL MPEXP (X, X)
      CALL MP40D (PLACES, X)
C COMPUTE EXP(PI/4)
      CALL MPDIVI (PI, 4, X)
      CALL MPEXP (X, X)
      CALL MP40D (PLACES, X)
C COMPUTE SIN(1) AND COS(1)
      CALL MPCIM (1, X)
      CALL MPCIS (X, X, Y, .TRUE.)
      CALL MP40D (PLACES, Y)
      CALL MP40D (PLACES, X)
C COMPUTE ZETA(3)
      CALL MPZETA (3, X)
      CALL MP40D (PLACES, X)
C COMPUTE LN(PHI), 1/LN(PHI), AND -LN(LN(2))
      CALL MPLN(PHI, X)
      CALL MP40D (PLACES, X)
      CALL MPREC (X, X)
      CALL MP40D (PLACES, X)
      CALL MPLNI (2, X)
      CALL MPLN (X, X)
      X(1) = -X(1)
      CALL MP40D (PLACES, X)
      WRITE (LUN, 80) MXSPTR
   80 FORMAT (/ 18H END OF TEST, USED, I5,
     $          23H WORDS OF WORKING SPACE //)
      STOP
      END
