c $$                   ******  test  ******
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 40 decimal places, but
c to increase the accuracy it is only necessary to change
c the statement places = 40, 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
c correct output (excluding headings) is as follows
c
c      1.4142135623 7309504880 1688724209 6980785697
c      1.7320508075 6887729352 7446341505 8723669428
c      2.2360679774 9978969640 9173668731 2762354406
c      3.1622776601 6837933199 8893544432 7185337196
c      1.2599210498 9487316476 7210607278 2283505703
c      1.4422495703 0740838232 1638310780 1095883919
c      1.1892071150 0272106671 7499970560 4759152930
c      0.6931471805 5994530941 7232121458 1765680755
c      1.0986122886 6810969139 5245236922 5257046475
c      2.3025850929 9404568401 7991454684 3642076011
c      1.4426950408 8896340735 9924681001 8921374266
c      0.4342944819 0325182765 1128918916 6050822944
c      3.1415926535 8979323846 2643383279 5028841972
c      0.0174532925 1994329576 9236907684 8861271344
c      0.3183098861 8379067153 7767526745 0287240689
c      9.8696044010 8935861883 4490999876 1511353137
c      1.7724538509 0551602729 8167483341 1451827975
c      2.6789385347 0774763365 5692940974 6776441287
c      1.3541179394 2640041694 5288028154 5137855193
c      2.7182818284 5904523536 0287471352 6624977572
c      0.3678794411 7144232159 5523770161 4608674458
c      7.3890560989 3065022723 0427460575 0078131803
c      0.5772156649 0153286060 6512090082 4024310422
c      1.1447298858 4940017414 3427351353 0587116473
c      1.6180339887 4989484820 4586834365 6381177203
c      1.7810724179 9019798523 6504103107 1795491696
c      2.1932800507 3801545655 9769659278 7382234616
c      0.8414709848 0789650665 2502321630 2989996226
c      0.5403023058 6813971740 0936607442 9766037323
c      1.2020569031 5959428539 9738161511 4499907650
c      0.4812118250 5960344749 7758913424 3684231352
c      2.0780869212 3502753760 1322606117 7957677422
c      0.3665129205 8166432701 2439158232 6694694543
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(500)
c temporary mp variables require space t+2 and t .le. 25.
      integer i, places, x(27), y(27), phi(27), pi(27)
c
c set output unit = 6 and working precision to the
c equivalent of at least 43 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 = 40
      call mpset2 (6, places+3, 27, 1, 500)
      write (lun, 5) b, t
    5 format (29h1test of mp package,   base =, i9,
     $        11h,  digits =, i3 ///)
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, i4,
     $          23h words of working space //)
      stop
      end