c $$                   ******  test2  ******
c test2 program for mp package
c
c this program tests various mp routines, especially those not
c called by program test.  it computes the constants given in
c computer approximations (by hart, cheney, lawson, maehly,
c mesztenyi, rice, thacher and witzgall, john wiley, 1968),
c appendix c, pp. 182-183, and various other constants
c which are described in the comments below.  the constants
c are computed to 40 significant figures, with working precision
c equivalent to at least 43 significant figures.
c
c correct output on a univac 1100/42 (36-bit word) is
c as follows.  on machines with wordlength other than 36 bits
c there will be some minor differences.  the results given
c below are correctly rounded to 40 significant figures.
c
c test2 of mp package,   base =    65536,  digits = 10
c
c constants in hart et al (order different)
c
c  4.848136811095359935899141023579479759564e-0005
c  1.745329251994329576923690768488612713443e-0002
c  3.926990816987241548078304229099378605246e-0001
c  5.641895835477562869480794515607725858441e-0001
c  6.366197723675813430755350534900574481378e-0001
c  7.853981633974483096156608458198757210493e-0001
c  7.978845608028653558798921198687637369517e-0001
c  9.189385332046727417803297364056176398614e-0001
c  1.570796326794896619231321691639751442099e+0000
c  1.772453850905516027298167483341145182798e+0000
c  2.356194490192344928846982537459627163148e+0000
c  3.141592653589793238462643383279502884197e+0000
c  6.283185307179586476925286766559005768394e+0000
c  3.678794411714423215955237701614608674458e-0001
c  7.788007830714048682451702669783206472968e-0001
c  1.284025416687741484073420568062436458336e+0000
c  2.718281828459045235360287471352662497757e+0000
c  5.772156649015328606065120900824024310422e-0001
c  1.090507732665257659207010655760707978993e+0000
c  1.189207115002721066717499970560475915293e+0000
c  1.414213562373095048801688724209698078570e+0000
c  3.162277660168379331998893544432718533720e+0000
c  1.259921049894873164767210607278228350570e+0000
c  1.587401051968199474751705639272308260391e+0000
c  2.154434690031883721759293566519350495259e+0000
c  4.641588833612778892410076350919446576551e+0000
c  1.442695040888963407359924681001892137427e+0000
c  3.321928094887362347870319429489390175865e+0000
c  3.465735902799726547086160607290882840378e-0001
c  6.931471805599453094172321214581765680755e-0001
c  1.386294361119890618834464242916353136151e+0000
c  2.302585092994045684017991454684364207601e+0000
c  3.010299956639811952137388947244930267682e-0001
c  4.342944819032518276511289189166050822944e-0001
c  1.305261922200515915484062278954890101937e-0001
c  1.950903220161282678482848684770222409277e-0001
c  2.588190451025207623488988376240483283491e-0001
c  5.000000000000000000000000000000000000000e-0001
c  7.071067811865475244008443621048490392848e-0001
c  8.660254037844386467637231707529361834714e-0001
c  2.474039592545229295968487048493891958934e-0001
c  4.794255386042030002732879352155713880818e-0001
c  8.414709848078965066525023216302989996226e-0001
c  9.238795325112867561281831893967882868224e-0001
c  9.659258262890682867497431997288973676339e-0001
c  9.849140335716425307719752129132743229305e-0002
c  1.989123673796580069115976226446762285979e-0001
c  2.679491924311227064725536584941276330572e-0001
c  4.142135623730950488016887242096980785697e-0001
c  5.773502691896257645091487805019574556476e-0001
c  1.000000000000000000000000000000000000000e+0000
c  1.732050807568877293527446341505872366943e+0000
c  2.553419212210362665044822364904736782042e-0001
c  5.463024898437905132551794657802853832976e-0001
c  1.557407724654902230506974807458360173087e+0000
c  2.864256869924070484773772302378002636590e-0001
c  4.563200724962071954298288879197037326544e-0001
c  5.323280990604765985693418601514840681969e-0001
c  6.020994932348296742913495709025331582540e-0001
c  7.236648875473587964064629355277428962481e-0001
c  7.761882222734748449391348929091469265099e-0001
c  8.238049667103636277653479963181034453800e-0001
c  8.669840139434846676407863266118724492986e-0001
c  9.061778032186181584326793796610085912753e-0001
c  9.184307960450497544469833659648160928434e-0001
c  9.418086078468580771196106632361044035158e-0001
c  9.742619823209360312818256760119330624856e-0001
c  8.516319137048080127004060150609260682003e-0001
c  4.720012157682347674476683878725009623642e-0001
c  3.631878383468673317955937477889247216476e-0001
c  5.668240889058739377112449634671602835403e-0001
c  2.369974904082422018721147551606796861398e-0011
c  2.326614794865976450546482207237974647586e-0008
c  4.769362762044698733814183536431305598090e-0001
c
c test of mpasin and mpatan
c
c  1.000016667416711312562227707199038367857e-0002
c -5.235987755982988730771072305465838140329e-0001
c  1.429256853470469400485532334664724427105e+0000
c  9.999666686665238206340116209279548561369e-0003
c -6.435011087932843868028092287173226380415e-0001
c  1.471127674303734591852875571761730851855e+0000
c
c test of mpbesj
c
c  9.999750001562495659729003899468320681723e-0001
c  9.975015620660400322812868984747920848320e-0001
c  7.651976865579665514497175261026632209093e-0001
c -2.459357644513483351977608624853287538296e-0001
c  1.998585030422312242422839095084899068063e-0002
c  2.478668615242017456133073111569370878617e-0002
c  4.999937500260416124132622612282082222967e-0003
c  4.993752603624199755633655243780648405856e-0002
c  4.400505857449335159596822037189149131274e-0001
c  4.347274616886143666974876802585928830627e-0002
c -7.714535201411215803268549492723447021161e-0002
c  4.728311907089523917576071901216916285418e-0003
c  2.170131138404967281693651142150815094613e-0017
c  2.169363960376002380635265343042715360913e-0011
c  2.093833800238926996560701453800780000026e-0005
c -1.445884208478510531774561260148174874671e-0002
c -3.352538314417667427285301848429116213846e-0002
c -2.469801093420243996535541028052587615959e-0002
c  2.368598385174410068779331157930426751244e-0274
c  2.368519166467611069896185607645793014427e-0201
c  2.360610483197187577029325106563488361132e-0128
c  1.688254978074890509296827170625271452023e-0055
c  9.633817342050361437615884024435958149431e-0002
c  7.173642352660120247867318411619579007593e-0003
c  1.300924406738082128445548553616298220657e-0671
c  1.300904893017434872887808473358528762111e-0507
c  1.298954989483319596748137632393071785666e-0343
c  1.117943529279165433085063648308631625217e-0179
c  1.547555253553207391363466511764830851880e-0022
c  1.152259751253901478112859948236454927804e-0002
c  1.944998262170446054748125826967690995132e-3818
c  1.944992252353628101811176793079875375991e-3018
c  1.944391364322010861731787302600779281025e-2218
c  1.885229418999920133044032342021032255068e-1418
c  8.526917009536047015776970746427276966500e-0620
c -2.987223375566258955506328177917819328443e-0002
c
c test of mperf, mperfc and mpdaw
c
c -1.000000000000000000000000000000000000000e+0000
c  2.000000000000000000000000000000000000000e+0000
c -5.000250037509378282727375137642333908163e-0003
c -1.000000000000000000000000000000000000000e+0000
c  2.000000000000000000000000000000000000000e+0000
c -5.025384718759852803274841986071548588791e-0002
c -8.427007929497148693412206350826092592961e-0001
c  1.842700792949714869341220635082609259296e+0000
c -5.380795069127684191363874204075567547920e-0001
c -1.124629160182848922032750717439683832217e-0001
c  1.112462916018284892203275071743968383222e+0000
c -9.933599239785286114978869519231223541097e-0002
c  1.124629160182848922032750717439683832217e-0001
c  8.875370839817151077967249282560316167783e-0001
c  9.933599239785286114978869519231223541097e-0002
c  8.427007929497148693412206350826092592961e-0001
c  1.572992070502851306587793649173907407039e-0001
c  5.380795069127684191363874204075567547920e-0001
c  1.000000000000000000000000000000000000000e+0000
c  2.088487583762544757000786294957788611561e-0045
c  5.025384718759852803274841986071548588791e-0002
c  1.000000000000000000000000000000000000000e+0000
c  6.405961424921732039021339148586394148214e-4346
c  5.000250037509378282727375137642333908163e-0003
c
c test of mpgam
c
c  3.993599489724414140919056205879117156950e-0038
c  4.113526718873633300615263770795138677800e+0000
c  9.999422883231624190805737422564434215028e+0003
c  3.447019240352198953918716891440225225102e+0002
c  1.216451004088320000000000000000000000000e+0017
c  4.023872600770937735437024339230039857194e+2564
c  1.272301195695055464182244180377444569507e+2566
c
c test of mpsin, cos, tan, sinh, cosh, tanh
c
c  5.063656411097587936565576104597854320650e-0001
c  8.623188722876839341019385139508425355101e-0001
c  5.872139151569290766778096356445878942588e-0001
c -1.344058570908067724206312775790006793681e+0043
c  1.344058570908067724206312775790006793681e+0043
c -1.000000000000000000000000000000000000000e+0000
c  5.440211108893698134047476618513772816836e-0001
c -8.390715290764524522588639478240648345199e-0001
c -6.483608274590866712591249330098086768169e-0001
c -1.101323287470339337723652455484636440290e+0004
c  1.101323292010332313972137609043787996345e+0004
c -9.999999958776927636195928371382757410508e-0001
c -8.414709848078965066525023216302989996226e-0001
c  5.403023058681397174009366074429766037323e-0001
c -1.557407724654902230506974807458360173087e+0000
c -1.175201193643801456882381850595600815156e+0000
c  1.543080634815243778477905620757061682602e+0000
c -7.615941559557648881194582826047935904128e-0001
c -9.983341664682815230681419841062202698992e-0002
c  9.950041652780257660955619878038702948386e-0001
c -1.003346720854505450580800457811115368190e-0001
c -1.001667500198440258237293835219050235149e-0001
c  1.005004168055803598987978442968341644710e+0000
c -9.966799462495581711830508367835218353896e-0002
c  9.983341664682815230681419841062202698992e-0002
c  9.950041652780257660955619878038702948386e-0001
c  1.003346720854505450580800457811115368190e-0001
c  1.001667500198440258237293835219050235149e-0001
c  1.005004168055803598987978442968341644710e+0000
c  9.966799462495581711830508367835218353896e-0002
c  8.414709848078965066525023216302989996226e-0001
c  5.403023058681397174009366074429766037323e-0001
c  1.557407724654902230506974807458360173087e+0000
c  1.175201193643801456882381850595600815156e+0000
c  1.543080634815243778477905620757061682602e+0000
c  7.615941559557648881194582826047935904128e-0001
c -5.440211108893698134047476618513772816836e-0001
c -8.390715290764524522588639478240648345199e-0001
c  6.483608274590866712591249330098086768169e-0001
c  1.101323287470339337723652455484636440290e+0004
c  1.101323292010332313972137609043787996345e+0004
c  9.999999958776927636195928371382757410508e-0001
c -5.063656411097587936565576104597854320650e-0001
c  8.623188722876839341019385139508425355101e-0001
c -5.872139151569290766778096356445878942588e-0001
c  1.344058570908067724206312775790006793681e+0043
c  1.344058570908067724206312775790006793681e+0043
c  1.000000000000000000000000000000000000000e+0000
c
c test of mpli and mpei
c
c -1.247678559626375088840194710872203816285e-0004
c -6.353279339727591513585474237270429058630e+0000
c -6.355232464831071802614455519358032212938e+0000
c -6.354744461697291123056164840748548577634e+0000
c  1.892463243835693144140481263216406702938e+0000
c -2.197435427851990963896917162334182400288e-0001
c -6.353767899171421044877727316545766253157e+0000
c  1.897772388875768322684241430851491375948e+0000
c -2.190250276806367219674691602503337494018e-0001
c -5.299113807150346461381605153457779262963e-0004
c  3.156552081410363613766570878788022081110e+0000
c -1.075880501319490268122213298905290428959e-0001
c  6.165599504787297937522981752669522749131e+0000
c  2.492228976241877759138440143998524848990e+0003
c -4.156968929685324277402859810278180384346e-0006
c  3.012614158407962992590174133903218497960e+0001
c  2.715552744853879821914014642310825410296e+0041
c -3.683597761682032180235192620508118987655e-0046
c  1.776096579901522266876406239486993179786e+0002
c  1.972045137141238302809645048412023552690e+0431
c -5.070893060235166549927200999685925144667e-0438
c  1.246137215899388459692771107529059792487e+0003
c  8.807699083674714448970900245119101084640e+4338
c -1.135370339631071751833431428951800778261e-4347
c
c test of mpzeta
c
c  1.202056903159594285399738161511449990765e+0000
c  1.082323233711138191516003696541167902775e+0000
c  1.036927755143369926331365486457034168057e+0000
c  1.000994575127818085337145958900319017006e+0000
c  1.000000953962033872796113152038683449346e+0000
c  1.000000000000909494784026388928253311839e+0000
c
c test of mpcam, mpgcda and mpgcdb
c
c  1.230000000000000000000000000000000000000e-0004
c -5.000000000000000000000000000000000000000e+1232
c  2.345600000000000000000000000000000000000e+0790
c  2.ef405befc1ce7a4cd25b45dc41336$+00000290
c  1.524157875319616034331961603431672002468e+39
c -2.895899851996160343319616034316720024680e+39
c  1.234567890120000000000000000000000000000e+11
c  1.234567890123456789012345678900000000000e+28
c -2.345678901234567890123456789000000000000e+28
c
c test varying digits, base and rounding rule
c
c  2.625374123e+17
c  2.625374124e+17
c  2.625374121e+17
c  2.625374132e+17
c  2.625374128e+17
c  2.625374128e+17
c  2.625374117e+17
c  2.625374129e+17
c  2.625374129e+17
c  2.625374130e+17
c  2.625374121e+17
c  2.625374133e+17
c  2.6253741264076874396e+17
c  2.6253741264076874396e+17
c  2.6253741264076874390e+17
c  2.6253741264076874416e+17
c  2.6253741264076874407e+17
c  2.6253741264076874409e+17
c  2.6253741264076874384e+17
c  2.6253741264076874410e+17
c  2.6253741264076874398e+17
c  2.6253741264076874400e+17
c  2.6253741264076874393e+17
c  2.6253741264076874419e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768744000000000000e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768744000000000000e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768743999999999999e+17
c  2.62537412640768744000000000000e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725971e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725971e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725972e+17
c  2.625374126407687439999999999992500725971e+17
c  2.625374126407687439999999999992500725972e+17
c
c test of various mp routines
c
c  1.144729885849400174143427351353058711647e+00
c  2.245915771836104547342715220454373502759e+01
c  3.095891371701158559100259628453669451579e+01
c -3.141592653589793238462643383279502884197e+00
c -3.141592653589793238462643383279502884197e+**
c -3.141592653589793238462643383279502884197e+1234
c  3.100627660751343     1
c  3.100627668029982014200000000000     1
c  3.1006276680299820148e+0001
c  5.915477618586312249935148194956040327282e+0039
c  1                     -8             1       655     23592     62914
c                                   36700     10485     49807     23592
c                                   62914     36701
c  1            -8589934590             1         0         0         0
c                                       0         0         0         0
c                                       0         0
c  1             8589934591         65535     65535     65535     65535
c                                   65535     65535     65535     65535
c                                   65535     65535
c        ************* ********** ********** **********
c  0.000000000000000000000000000000000000000e+0000
c
c final contents of common /mpcom/
c            65536              10      8589934591               6
c              500               1             266               1
c       8589934591     -8589934591               0               1
c                0              43              27     34359738367
c                6              80              10              10
c                e               6               1
c
c end of test2, used  266 words of working space
c
c to conform to the fortran standard we should declare blank and
c labelled common here, for example -
c
      common mpwork
      common /mpcom/ mppar
      integer mppar(23), mpwork(500)
c
c this is necessary on some systems (e.g. burroughs b6700),
c but not on most.  see comments in mpinit for more details.
c
c t .le. 25 for wordlength at least 16 bits.  dimensions of
c w, x etc. can be reduced if wordlength exceeds 16 bits.
c dimension of w, x etc. must agree with mt2 - see below.
c dimension of a must agree with idima - see below.
c
      integer a(43), base, decpl, i, id, idima, idummy(1), ilim,
     $  j, jl, j1(3), j10(3), j11(7), j12(7), j13(8), j14(4),
     $  j15(6), j16(6), j2(4), j3(3), j4(4), j5(6), j6(7), j7(12),
     $  j8(3), j9(3), k, ktu, ktv, lun, mpc(23), mt2, mxsptr, n,
     $  numdig, n2, pi(27), sv, w(27), x(27), y(27), z(27),
     $  mpchgb, mpdigs, mpparn, mpsiga
      logical error
      real rx
      double precision dx
c
c mt2 must equal dimension of arrays w, x etc. above (at least
c t+2 for t-digit mp numbers).
c
c idima must equal dimension of array a above (at least decpl+3
c for decpl significant decimal places in output).
c
      data mt2 /27/,  idima /43/
c
c data for test arguments
c
      data j1(1), j1(2), j1(3) /16200, 45, 2/
      data j2(1), j2(2), j2(3), j2(4) /-1, -4, 4, 1/
      data j3(1), j3(2), j3(3) /8, 4, 2/
      data j4(1), j4(2), j4(3), j4(4) /2, 4, 10, 100/
      data j5(1), j5(2), j5(3) /24, 16, 12/
      data j5(4), j5(5), j5(6) /6, 4, 3/
      data j6(1), j6(2), j6(3) /32, 16, 12/
      data j6(4), j6(5), j6(6), j6(7) /8, 6, 4, 3/
      data j7(1), j7(2), j7(3), j7(4) /9, 15, 18, 21/
      data j7(5), j7(6), j7(7), j7(8) /27, 30, 33, 36/
      data j7(9), j7(10), j7(11), j7(12) /39, 40, 42, 45/
      data j8(1), j8(2), j8(3) /0, 1, 10/
      data j9(1), j9(2), j9(3) /1, -50, 99/
      data j10(1), j10(2), j10(3) /4, -300, 4000/
      data j11(1), j11(2), j11(3), j11(4) /-101, -13, 1, 33/
      data j11(5), j11(6), j11(7) /20, 1000, 2001/
      data j12(1), j12(2), j12(3) /3, 7, 10000/
      data j12(4), j12(5), j12(6), j12(7) /5, 1, 1, 2/
      data j13(1), j13(2), j13(3) /1, 1023, 1025/
      data j13(4), j13(5), j13(6) /1486, 10, 100/
      data j13(7), j13(8) /1000, 10000/
      data j14(1), j14(2), j14(3), j14(4) /-1, 3, -4, 2/
      data j15(1), j15(2), j15(3) /0, 1, 6/
      data j15(4), j15(5), j15(6) /73, 164, 800/
      data j16(1), j16(2), j16(3) /3, 4, 5/
      data j16(4), j16(5), j16(6) /10, 20, 40/
c
c initialize using mpinit.  the standard version of mpinit sets
c the precision to the equivalent of at least 43 decimal places.
c if nonstandard versions of mpinit are used, the dimension and
c data statements above may have to be changed.
c
      call mpinit (idummy)
c
c use mpparm to get the base, number of digits etc.
c it is easier to use the function mppara, but the pfort
c verifier objects to functions with hollerith arguments.
c
      call mpparm (4hbase, .false., base)
      call mpparm (6hnumdig, .false., numdig)
      call mpparm (3hlun, .false., lun)
c
c see above regarding idima.
      call mpparm (5hdecpl, .false., decpl)
      decpl = min0 (idima - 3, decpl)
c
      write (lun, 10) base, numdig
   10 format (30h1test2 of mp package,   base =, i9,
     $  11h,  digits =, i3)
c
c check that version of mpinit used is compatible with test2.
c
      n2 = numdig + 2
      if (mt2 .ge. n2) go to 20
      write (lun, 15) n2
   15 format (48h *** increase mt2 and dimensions of w, x etc. to /
     $  9h at least, i6, 13h in test2 ***)
      call mperr
      stop
c
c compute pi using gauss-lagrange method
   20 call mppigl (pi)
c
c compute constants given in hart et al, in approximately the
c same order as given there.
      write (lun, 25)
   25 format (/ 42h constants in hart et al (order different) /)
c
c compute pi/64800, pi/180, pi/8
      call mpdivi (pi, 4, y)
      do 30 i = 1, 3
      call mpdivi (y, j1(i), x)
c write on unit lun (see call to mpset above).
   30 call mpfout (x, decpl)
c compute sqrt(1/pi)
      call mproot (pi, -2, x)
      call mpfout (x, decpl)
c compute 2/pi
      call mpdivi (pi, 2, z)
      call mprec (z, z)
      call mpfout (z, decpl)
c print pi/4 (already in y)
      call mpfout (y, decpl)
c compute sqrt (2/pi)
      call mpsqrt (z, x)
      call mpfout (x, decpl)
c compute ln(sqrt(2*pi))
      call mpmuli (pi, 2, z)
      call mpln (z, x)
      call mpdivi (x, 2, x)
      call mpfout (x, decpl)
c compute pi/2
      call mpdivi (pi, 2, x)
      call mpfout (x, decpl)
c compute sqrt (pi)
      call mpsqrt (pi, x)
      call mpfout (x, decpl)
c compute 3pi/4, pi, 2pi
      call mpmulq (pi, 3, 4, x)
      call mpfout (x, decpl)
      call mpfout (pi, decpl)
      call mpfout (z, decpl)
c compute exp(-1), exp(-1/4), exp(1/4), exp(1)
      do 35 i = 1, 4
      call mpcqm (1, j2(i), x)
      call mpexp (x, x)
   35 call mpfout (x, decpl)
c compute eulers constant
      call mpeul (x)
      call mpfout (x, decpl)
c compute sqrt(sqrt(sqrt(2))), sqrt(sqrt(2)), sqrt(2)
      do 40 i = 1, 3
      call mpqpwr (2, 1, 1, j3(i), x)
   40 call mpfout (x, decpl)
c compute sqrt(10)
      call mpqpwr (10, 1, 1, 2, x)
      call mpfout (x, decpl)
c compute cube root of 2, 4, 10, 100
      do 45 i = 1, 4
      call mpqpwr (j4(i), 1, 1, 3, x)
   45 call mpfout (x, decpl)
c compute log2(e), log2(10)
      call mplni (2, w)
      call mprec (w, y)
      call mpfout (y, decpl)
      call mplni (10, z)
      call mpmul (z, y, x)
      call mpfout (x, decpl)
c compute ln(sqrt(2)), ln(2), ln(4), ln(10)
      call mpdivi (w, 2, x)
      call mpfout (x, decpl)
      call mpfout (w, decpl)
      call mpmuli (w, 2, x)
      call mpfout (x, decpl)
      call mpfout (z, decpl)
c compute log10(2), log10(e)
      call mpcim (2, x)
      call mplg10 (x, x)
      call mpfout (x, decpl)
      call mprec (z, x)
      call mpfout (x, decpl)
c compute sin(pi/j) for j = 24, 16, 12, 6, 4, 3
c note - order is slightly different from hart et al here
      do 50 i = 1, 6
      call mpdivi (pi, j5(i), x)
      call mpsin (x, x)
   50 call mpfout (x, decpl)
c compute sin(1/4), sin(1/2), sin(1)
      do 60 i = 1, 3
      call mpcqm (1, 2**(3-i), x)
      call mpsin (x, x)
   60 call mpfout (x, decpl)
c compute sin(3pi/8), sin(5pi/12)
      call mpmulq (pi, 3, 8, x)
      call mpsin (x, x)
      call mpfout (x, decpl)
      call mpmulq (pi, 5, 12, x)
      call mpsin (x, x)
      call mpfout (x, decpl)
c compute tan (pi/j) for j = 32, 16, 12, 8, 6, 4, 3
c note - order is slightly different from hart et al
      do 70 i = 1, 7
      call mpdivi (pi, j6(i), x)
      call mptan (x, x)
   70 call mpfout (x, decpl)
c compute tan (1/4), tan (1/2), tan(1)
      do 80 i = 1, 3
      call mpcqm (1, 2**(3-i), x)
      call mptan (x, x)
   80 call mpfout (x, decpl)
c compute tan(j*pi/96) for j = 9, 15, 18, 21, 27, 30,
c                              33, 36, 39, 40, 42, 45
      call mpcim (96, y)
      do 90 i = 1, 12
      call mpmuli (pi, j7(i), x)
      call mpatn2 (x, y, x)
   90 call mpfout (x, decpl)
c compute j(nu, pi/4) and j(nu, pi/2) for nu = 0, 1, 10
      do 100 i = 1, 3
      call mpdivi (pi, 4, x)
      call mpbesj (x, j8(i), x)
      call mpfout (x, decpl)
      call mpdivi (pi, 2, x)
      call mpbesj (x, j8(i), x)
  100 call mpfout (x, decpl)
c compute arcerf(1/2) using newtons method with first
c approximation 0.477
      call mpsqrt (pi, y)
      call mpdivi (y, 2, y)
c compute number of iterations necessary
      ilim = mpchgb (2, (decpl+1)/2, 1)
      call mpcqm (477, 1000, x)
c could save time by reducing t at first (see eg mprec)
      do 110 i = 1, ilim
      call mperf (x, z)
      call mpaddq (z, -1, 2, z)
      call mpmul (x, x, w)
      call mpexp (w, w)
      call mpmul (w, y, w)
      call mpmul (w, z, z)
  110 call mpsub (x, z, x)
      call mpfout (x, decpl)
c
c finished with constants in hart et al, now compute some more
c to test other mp routines
c
      write (lun, 115)
  115 format (/ 26h test of mpasin and mpatan /)
c compute asin(1/100), asin(-1/2), asin(99/100)
      do 120 i = 1, 3
      call mpcqm (j9(i), 100, x)
      call mpasin (x, x)
  120 call mpfout (x, decpl)
c compute atan(1/100), atan(-3/4), atan(10)
      do 130 i = 1, 3
      call mpcqm (j10(i), 400, x)
      call mpatan (x, x)
  130 call mpfout (x, decpl)
      write (lun, 135)
  135 format (/ 15h test of mpbesj /)
c compute j(nu, x) for x = 0.01, 0.1, 1, 10, 100, 1000
c                  and nu = 0, 1, 6, 73, 164, 800
      do 150 i = 1, 6
      do 140 j = 1, 6
      if (j.lt.3) call mpcqm (1, 10**(3-j), x)
      if (j.ge.3) call mpcqm (10**(j-3), 1, x)
      call mpbesj (x, j15(i), x)
  140 call mpfout (x, decpl)
  150 continue
      write (lun, 155)
  155 format (/ 32h test of mperf, mperfc and mpdaw /)
c compute erf(x), erfc(x), and daw(x) (dawsons integral) for
c x = -100, -10, -1, -0.1, 0.1, 1, 10, 100
      do 160 i = 1, 8
      if (i.le.4) call mpcqm (10**(4-i), -10, x)
      if (i.gt.4) call mpcqm (10**(i-5), 10, x)
      call mperf (x, y)
      call mpfout (y, decpl)
      call mperfc (x, y)
      call mpfout (y, decpl)
      call mpdaw (x, y)
  160 call mpfout (y, decpl)
      write (lun, 165)
  165 format (/ 14h test of mpgam /)
c compute gamma(x) for x = -101/3, -13/7, 1/10000, 33/5, 20,
c                          1000, 2001/2
      do 170 i = 1, 7
      call mpcqm (j11(i), j12(i), x)
      call mpgam (x, x)
  170 call mpfout (x, decpl)
      write (lun, 175)
  175 format (/ 42h test of mpsin, cos, tan, sinh, cosh, tanh /)
c compute sin(x), cos(x), tan(x), sinh(x), cosh(x) and tanh(x)
c for x = -100, -10, -1, -0.1, 0.1, 1, 10, 100
      do 180 i = 1, 8
      if (i.le.4) call mpcqm (10**(4-i), -10, x)
      if (i.gt.4) call mpcqm (10**(i-5), 10, x)
      call mpsin (x, y)
      call mpfout (y, decpl)
      call mpcos (x, y)
      call mpfout (y, decpl)
      call mptan (x, y)
      call mpfout (y, decpl)
      call mpsinh (x, y)
      call mpfout (y, decpl)
      call mpcosh (x, y)
      call mpfout (y, decpl)
      call mptanh (x, y)
  180 call mpfout (y, decpl)
      write (lun, 185)
  185 format (/ 22h test of mpli and mpei /)
c compute li(x), ei(x), ei(-x) for x = 1/1024, 1023/1024,
c 1025/1024, 1486/1024, 10, 100, 1000, 10000
      do 190 i = 1, 8
      id = 1024
      if (i.gt.4) id = 1
      call mpcqm (j13(i), id, x)
      call mpli (x, y)
      call mpfout (y, decpl)
      call mpei (x, y)
      call mpfout (y, decpl)
      call mpneg (x, x)
      call mpei (x, y)
  190 call mpfout (y, decpl)
      write (lun, 195)
  195 format (/ 15h test of mpzeta /)
c compute zeta(i) for i = 3, 4, 5, 10, 20, 40
      do 200 i = 1, 6
      call mpzeta (j16(i), x)
  200 call mpfout (x, decpl)
      write (lun, 210)
  210 format (/ 33h test of mpcam, mpgcda and mpgcdb /)
c test mpcam.
      call mpcam (7h123e-6$, y)
      call mpfout (y, decpl)
      call mpcam (10h-.05d1234$, y)
      call mpfout (y, decpl)
      call mpcam (12h+23.456+789$, y)
      call mpfout (y, decpl)
c change output base to 16, exponent field width to 10
      call mpparb (16, 7houtbase)
      call mpparb (10, 7hexwidth)
c write y again, this time with base 16 and 30 digits
      call mpfout (y, 30)
c restore output base to 10, set exponent width to 4
      call mpparb (10, 3hout)
      call mpparb (4, 3hexw)
      call mpcam (43h1524157875319616034331961603431672002468.0$, x)
      call mpfout (x, decpl)
      call mpcam (42h-2895899851996160343319616034316720024680$, y)
      call mpfout (y, decpl)
c test mpgcda
      call mpgcda (x, y, z)
      call mpfout (z, decpl)
c test mpgcdb
      call mpgcdb (x, y)
      call mpfout (x, decpl)
      call mpfout (y, decpl)
      write (lun, 220)
  220 format (/ 44h test varying digits, base and rounding rule /)
c save base, digits and rounding indicator, then vary them
      call mpsavn (sv)
      jl = max0 (2, base - 2)
      do 230 i = 1, 4
      do 230 j = jl, base
c set base to j
      call mpparb (j, 4hbase)
c set number of digits to equivalent of at least 10*i decimal digits
      call mpparb (mpdigs (10*i), 6hnumdig)
      do 230 k = 1, 4
c set rndrl to k-1
      call mpsetr (k-1)
c compute exp(pi*sqrt(163))
      call mpcam (2hpi, x)
      call mpqpwr (163, 1, 1, 2, y)
      call mpmul (x, y, x)
      call mpexp (x, x)
  230 call mpfout (x, 10*i)
c restore base, number of digits, and rndrl.
      call mpresn (sv)
      write (lun, 240)
  240 format (/ 28h test of various mp routines /)
c compute ln(pi) using mplngs
      call mplngs (pi, x)
      call mpfout (x, decpl)
c compute pi**e using mppwr2
      call mpcim (1, x)
      call mpexp (x, x)
      call mppwr2 (pi, x, x)
      call mpfout (x, decpl)
c compute 2*pi**3 - 4*pi**2 + 3*pi - 1 using mppoly
      call mppoly (pi, x, j14, 4)
      call mpfout (x, decpl)
c convert -pi to character format using mpout, then back to mp
c format using mpin.
      call mpneg (pi, x)
      call mpout (x, a, decpl+3, decpl)
      call mpin (a, x, decpl+3, error)
c error should be false
      if (error) call mperr
      call mpfout (x, decpl)
c x should be negative (if not, probable cause is minus and
c blank not distinguished by function mpis).
      if (mpsiga (x) .ge. 0) call mperr
c now test mpine by forming -pi*(10**1234)
      call mpine (a, x, decpl+3, 1234, error)
      if (error) call mperr
c exponent will print as e+**
      call mpfout (x, decpl)
c try again with larger exponent field
      call mpparb (6, 5hexwid)
      call mpfout (x, decpl)
c convert pi**3 to exponent and single-precision fraction
      call mppwr (pi, 3, x)
      call mpcmre (x, n, rx)
      write (lun, 250) rx, n
  250 format (1x, f18.15, i6)
c now convert to exponent and double-precision fraction
      call mpcmde (x, n, dx)
      write (lun, 260) dx, n
  260 format (1x, f33.30, i6)
c convert back to multiple-precision
      call mpcdm (dx*(10d0**n), x)
c and print (number of places agreement depends on
c floating-point fraction length)
      call mpfout (x, 20)
c compute integer part of pi**80
      call mppwr (pi, 80, x)
      call mpcmim (x, x)
      call mpfout (x, decpl)
c dump machine-precision, minreal and maxreal (these depend
c on base and number of digits, and hence on the wordlength
c of the machine used).
      call mpcam (7hepsilon, x)
      call mpdump (x)
      call mpcam (7hminreal, x)
      call mpdump (x)
      call mpcam (7hmaxreal, y)
      call mpdump (y)
c print using mp40d (y is too large so mpout gives
c all asterisks)
      call mp40d (decpl, y)
c cause an mp underflow (result z will be set to zero by mpunfl)
      call mpparm (6hktunfl, .false., ktu)
      call mpmul (x, x, z)
      call mpfout (z, decpl)
      call mpparm (6hktunfl, .false., ktv)
      ktu = ktv - ktu
c ktu should be one as one underflow should have occurred.
      if (ktu.ne.1) call mperr
c move mp parameters to local array
      do 270 i = 1, 23
  270 mpc (i) = mpparn (i)
c and write them
      write (lun, 280) mpc
  280 format (/ 33h final contents of common /mpcom/,
     $  5(/1x, 4i16) / 16x, a1, 2i16)
c get maximum stack pointer
      call mpparm (6hmxsptr, .false., mxsptr)
      write (lun, 290) mxsptr
  290 format (/ 19h end of test2, used, i5,
     $  23h words of working space /)
      stop
      end