1MP USER'S GUIDE (FOURTH EDITION)
0THE CHARACTER SET USED FOLLOWS, WITH  SPECIAL  CHARACTERS  (QUOTE,  LESS
 THAN,  GREATER  THAN,  UNDERLINE)  AND LOWER CASE LETTERS UNDERNEATH THE
 CORRESPONDING UPPER CASE ONES.
0ABCDEFGHIJKLMNOPQRSTUVWXYZ$0123456789+-*/=(),. QUOTE LT GT UNDERLINE
 abcdefghijklmnopqrstuvwxyz                       '   <  >      _
0This  file  is  optionally  distributed  in  upper-case  only, using the
 standard Fortran character set plus QUOTE ('), LESS THAN (<) and GREATER
 THAN (>).
0The standard Fortran carriage-control characters are in column  1.  Only
 blank,  '1', '+' and '0' are used.  There are at most 57 lines per page,
 and at most 72 characters (excluding the carriage-control character  and
 trailing blanks) per line.
1
0
0
0
0
0
0                  MP USER'S GUIDE  (Fourth Edition)
+                  _________________________________
0
0
0                          Richard P. Brent
0
0
0                   Department of Computer Science
0                   Australian National University
0                     Box 4, Canberra, ACT 2600
0                             Australia
0
0
0                    Technical Report TR-CS-81-08
0                             June 1981
0
0
0                 Copyright <c> 1981 by R. P. Brent
1                             MP User's Guide                     Page  2
+                             _______________
0CONTENTS
+________
0
       1. General description of MP                               page  3
0           1.1 Introduction                                            3
            1.2 Restrictions and assumptions                            4
            1.3 Note on correctness                                     5
            1.4 History and likely future development                   5
            1.5 Summary of useful MP routines                           7
0
       2. EXAMPLE program                                               8
0
       3. Parameters in COMMON /MPCOM/                                 10
0           3.1 Description of the parameters                          10
            3.2 Methods of setting the parameters                      12
            3.3 Efficiency considerations                              13
0      4. The Augment interface                                        14
0           4.1 Introduction                                           14
            4.2 The description deck                                   15
            4.3 JACOBI - an example of the use of Augment              18
0      5. Summary for Augment users                                    20
0           5.1 Notes                                                  20
            5.2 Summary                                                21
            5.3 Synonyms                                               25
            5.4 Reserved words                                         25
0
       6. Description of MP routines                                   26
0
       7. Test programs                                                68
0           7.1 TEST                                                   68
            7.2 TEST2                                                  68
            7.3 Other test programs                                    68
0      8. Comparison with earlier versions of MP                       69
0           8.1 New routines and capabilities                          69
            8.2 Incompatibilities with earlier versions of MP          70
0      9. Installation instructions                                    71
0           9.1 Essentials                                             71
            9.2 Conversion notes                                       72
            9.3 Desirable changes                                      73
            9.4 Note on Fortran 77                                     73
1Sec. 1.1                     MP User's Guide                     Page  3
+                             _______________
01. GENERAL DESCRIPTION OF MP
+____________________________
01.1 Introduction
+________________
0MP is a multiple-precision floating-point  arithmetic  package.   It  is
 almost  completely  machine-independent,  and  should run on any machine
 with an ANSI  Standard  Fortran  (ANS  X3.9-1966)  compiler,  sufficient
 memory,  and  a wordlength (for integer arithmetic) of at least 16 bits.
 A precompiler (Augment) which facilitates the use of the MP  package  is
 available - see Section 4 for details.
0MP has been tested on various machines including a Univac 1108, a Univac
 1100/82,  a  DEC 10, an IBM 360/50, 360/91 and 370/168, and a PDP 11/45,
 with various Fortran compilers.  These machines have  effective  integer
 wordlengths  ranging  from 16 to 48 bits.  A wordlength of 12 bits would
 probably be sufficient, but this has not been tested.
0MP  works  with  normalized  floating-point  numbers.   The base (B) and
 number of digits (T) are arbitrary, subject to some  restrictions  given
 below, and may be varied dynamically.
0T-digit floating-point numbers are stored in integer arrays of dimension
 at least T+2, with the following conventions
0      word 1 = sign (0, -1 or +1),
       word 2 = exponent (to base B),
       words 3 to T+2 = normalized fraction  (one base B digit per word).
0Note  that  words  2  to T+2 are undefined if the sign is zero.  Integer
 arrays used for storage of T-digit floating-point numbers in  the  above
 format will be referred to as (unpacked) multiple-precision variables or
 (unpacked)  multiple-precision  numbers.   There  is   also   a   packed
 representation,  where  two  base B digits are stored in each word - for
 details see the description of MPPACK in Section 6.
0The  user  may ask for rounded or truncated arithmetic to be used.  With
 rounded  arithmetic   the   basic   arithmetic   operations   (addition,
 multiplication and division)  give the correctly rounded T-digit, base B
 result.  With truncated arithmetic the relative error in the  result  is
 at  most  B**(1-T).  Directed roundings (useful for interval arithmetic)
 are also available - for details see Section 3.
0Rounding  options  are  implemented  for  some  elementary  and  special
 functions (e.g. MPSQRT) but not yet for others.  For details see Section
 6.    Most routines give  at  worst  a  result
 y  =  f(x)  which  could  be  obtained by making an O(B**(1-T)) relative
 perturbation in the argument x, evaluating the function f exactly,  then
 making an O(B**(1-T)) relative perturbation in the result.
0Exponents can lie in the range 1-M, ... , +M inclusive, where M  is  set
 by   the   user.    Thus,  numbers  from  B**(-M)  to  almost  B**M  are
 representable.  On underflow during an arithmetic operation, the  result
 is  set  to  zero by subroutine MPUNFL.  (This is the default action and
 may be modified - see Section 6.)   On  overflow  subroutine  MPOVFL  is
 called and execution is terminated with an error message.
1Sec. 1.2                     MP User's Guide                     Page  4
+                             _______________
0Error messages are printed on logical unit LUN, where LUN is set by  the
 user, and then execution is terminated by a call to subroutine MPERR.
0The parameters B, T, M, LUN etc. are stored in a labelled COMMON  block,
 COMMON  /MPCOM/.  A working array of sufficient size must be provided in
 blank COMMON.  Details are given in Section 3.1.
01.2 Restrictions and assumptions
+________________________________
0B (the base) must be at least 2.
0T (the number of digits) must be at least 2.
0M  (the  exponent range) must be greater than T and less than 1/4 of the
 largest machine-representable integer.
08*B**2-1 must not exceed the largest machine-representable integer.
0B**(T-1) should be at least 10**7 (this restriction does  not  apply  to
 the basic arithmetic operations).
0B, T etc. must be initialized by calling either MPINIT or MPSET2  before
 calling any other MP routines.  For details see Section 6.
0Blank COMMON rather than labelled COMMON is  used  for  working  storage
 because  a  curious  restriction  in  the  ANSI  (1966) Fortran standard
 requires that a labelled COMMON block be declared with the  same  length
 in  each  subprogram  in  which  it  is  declared,  but this restriction
 fortunately does not apply to blank COMMON.   Although  space  in  blank
 COMMON  must  be provided for the MP routines to use as working storage,
 this space need not start at the first word of  blank  COMMON  (see  the
 description  of  MPSET2 in Section 6).  Thus, programs which already use
 blank COMMON can make use of the MP package.
0MP  variables  used  as  arguments  of subroutines need not be distinct.
 Thus
       CALL MPADD (X, X, X)
 and
       CALL MPEXP (X, X)
0are  allowable.   However,  distinct  arrays which overlap should not be
 used.  For example,  CALL MPEXP (X(1), X(2)) is not allowed if X  is  an
 integer array holding a multiple-precision variable.
0
 It is assumed that the compiler  passes  addresses  of  arrays  used  as
 arguments  in  subroutine  calls  (i.e. call by reference), and does not
 check for array bounds violations (either for arguments  of  subroutines
 or for arrays in COMMON).  Apart from these violations, MP is written in
 a portable subset of ANSI Standard Fortran (ANSI X3.9-1966).   This  has
 been  checked by the PFORT verifier.  The only machine-dependent routine
 is MPUPW (which unpacks characters stored several  to  a  word).   Other
 routines  which  may require trivial changes are MPINIT, MPIS and MPLARG
 (see comments in Sections 6 and 9 below).  For comments on  Fortran  77,
 see Section 9.4.
1Sec. 1.3                     MP User's Guide                     Page  5
+                             _______________
01.3 Note on correctness
+_______________________
0It is clearly impossible to thoroughly test even one MP routine with all
 possible  combinations  of  B,  T,  M  etc.  It is also impracticable to
 formally prove correctness of any nontrivial MP routines  using  present
 theorem-proving  techniques.  Thus, no warranty can be given that the MP
 package is correct.  To gain confidence in it we rely on
0      1. Careful modular design, using the idea of 'information  hiding'
          and  the  avoidance  of  'tricks'  and 'side-effects' as far as
          possible.
0      2. Error analysis (in the style of Wilkinson).
0      3. Limited  testing  including,  where  possible,  comparison   of
          results  obtained  by  at  least  two independent methods (e.g.
          MPLN, MPLNGS and MPLNI compute logarithms by three  independent
          methods).
0Although no absolute assurance is possible, a user may be confident that
 a result obtained using MP is correct if
0      1. He uses more digits than are strictly necessary.
0      2. He repeats the computation using a different base B.
0      3. He  uses  two  independent  mathematical  formulations  of  the
          problem.
0It  is  unlikely  that  a bug in MP, the user's Fortran compiler, or the
 user's hardware, would remain undetected if steps 1-3 were followed.
0MP  makes  minimal  assumptions  - primarily that INTEGER arithmetic (on
 numbers which are not too large) is exact.
0MP  does  not use Fortran library routines such as SIN, ALOG, EXP, SQRT,
 and uses only 'A1' format for output of multiple-precision numbers.
0REAL  and  DOUBLE  PRECISION  arithmetic  are  used only in routines for
 conversion of these types to and from multiple-precision etc.
0
 1.4 History and likely future development
+_________________________________________
0The first working version of MP (version 731101) was written in November
 1973.  Between 1973 and 1978 a considerable number of  improvements  and
 additions were made, and nonstandard Fortran constructs were eliminated.
 In April 1978 the Augment interface  routines  (see  Section  4.1)  were
 added and the first version of the Augment description deck was written.
 The version of MP published in TOMS and formerly available from the  ACM
 Algorithms Distribution service (version 770217) did not include Augment
 interface routines, but a version of MP  which  includes  them  (version
 780802)  is now available from the Algorithms Distribution service.  The
 most recent version of MP may be obtained from the author at the address
 given on the title page of this User's Guide.
1Sec. 1.4                     MP User's Guide                     Page  6
+                             _______________
0In 1979 the storage allocation scheme  was  improved,  rounding  options
 were  implemented  for  many  MP routines, dependence on Fortran library
 routines such as EXP and ALOG and on  REAL  arithmetic  was  eliminated,
 many  routines  were  modified  to  allow  packed  as  well  as unpacked
 multiple-precision arguments, and corresponding changes were made in the
 Augment  description  deck.   For more details, see Section 8.1.  The MP
 User's Guide (Second Edition) was extensively  revised  to  produce  the
 third  edition, and at the same time it was converted from upper-case to
 duo-case.  Only minor changes and corrections have been made since then.
0For  an  introduction  to  the  design and philosophy of MP, see - R. P.
 Brent, 'A Fortran multiple-precision  arithmetic  package',  ACM  Trans.
 Math.   Software 4 (March 1978), 57-70.  Additional details are given in
 'Algorithm 524 - MP, a Fortran multiple-precision  arithmetic  package',
 ibid, 71-81, and in Section 6 below.  The mathematical theory behind the
 algorithms used in the MP  package  is  described  in  -  R.  P.  Brent,
 'Unrestricted   algorithms   for   elementary  and  special  functions',
 Information Processing 80, North Holland, 1980, 613-619.
0As  mentioned  above, a preprocessor (Augment) which facilitates the use
 of the MP package is available.   See  -  R. P. Brent, J. A. Hooper  and
 J.  M.  Yohe,  'An  Augment  interface  for  Brent's  multiple-precision
 arithmetic package', ACM Trans. Math. Software 6  (1980),  146-149,  and
 Sections 4 and 5 below.
0Motivation for certain details of our implementation of T-digit, base  B
 arithmetic  is given in - C. H. Reinsch, 'Principles and preferences for
 computer arithmetic', ACM Signum Newsletter 14, 1 (March  1979),  12-27.
 In  designing  the Augment interface we have been influenced by the IFIP
 Working Group 2.5 proposals - see B. Ford  (ed.),  'Parameterization  of
 the  environment for transportable numerical software', ACM Trans. Math.
 Software 4 (1978), 100-103.  See also - W. S. Brown and S.  I.  Feldman,
 'Environment   parameters   and   basic   functions  for  floating-point
 computations', ACM Signum Newsletter 14, 1 (March 1979), 42-45.
0In  the  future  we  hope  to  implement  rounding  options  for more MP
 routines, and write a  multiple-precision  interval  arithmetic  package
 which  uses MP and takes advantage of the directed rounding options.  An
 Augment interface to facilitate  the  use  of  the  interval  arithmetic
 package  would  be  provided.   (A similar project has been attempted by
 John P. Jeter -  see  his  'A  variable  precision  interval  data  type
 extension  to  Fortran',  M.  Sc. thesis, Univ. of South West Louisiana,
 July 1979.)
0A  never-ending  project  is to implement multiple-precision versions of
 ever more special functions, and to  improve  the  efficiency  of  those
 multiple-precision routines already implemented.
0Correspondence concerning MP should be addressed to Prof. R. P. Brent at
 the  address  given  on the title page.  Reports of bugs, information on
 applications for which the MP package has proved useful, and suggestions
 for the future development of the package are particularly welcome.
1Sec. 1.5                     MP User's Guide                     Page  7
+                             _______________
01.5 Summary of useful MP routines
+_________________________________
0We mention here the names of those MP routines which are likely to be of
 interest  to someone using the MP package without the aid of the Augment
 interface.  Routines which are called  by  other  MP  routines  but  are
 unlikely  to  be  called  directly by a user of the package are omitted.
 For the user of Augment, the summary given in Section 5.2 will  be  more
 useful.  Details of all MP routines may be found in Section 6.
0Basic arithmetic
0      MPADD, MPADDI,  MPADDQ,  MPDIV,  MPDIVI,  MPMUL,  MPMULI,  MPMULQ,
       MPREC, MPSUB
0Powers and roots
0      MPPWR, MPPWR2, MPQPWR, MPROOT, MPSQRT
0Elementary functions
0      MPASIN,  MPATAN,  MPATN2,  MPCIS,  MPCOS,  MPCOSH,  MPEXP, MPEXP1,
       MPLG10, MPLN, MPLNI, MPLNS, MPSIN, MPSINH, MPTAN, MPTANH
0Special functions
0      MPBESJ, MPDAW, MPEI, MPERF, MPERFC, MPGAM, MPGAMQ, MPLI, MPLNGM
0Constants
0      MPBERN, MPEPS, MPEUL, MPMAXR, MPMINR, MPPI, MPZETA
0Input and output
0      MPFIN, MPFOUT, MPIN, MPOUT, MPUNFR, MPUNFW
0Conversion
0      MPCAM, MPCDM, MPCIM, MPCMD, MPCMDE, MPCMI, MPCMR,  MPCMRE,  MPCQM,
       MPCRM
0Comparison
0      MPCMPA,  MPCMPD, MPCMPI, MPCMPQ, MPCMPR, MPCOMP, MPEQ, MPGE, MPGT,
       MPLE, MPLT, MPNE
0General utility routines
0      MPABS, MPCEIL,  MPCHEB,  MPCHEV,  MPCMF,  MPCMIM,  MPDIGS,  MPDIM,
       MPFLOR,  MPGCDA,  MPGCDB,  MPINIT,  MPMAX,  MPMIN,  MPMOD,  MPNEG,
       MPPACK, MPPARA, MPPARB, MPPOLY,  MPSETR,  MPSET2,  MPSIGN,  MPSTR,
       MPUNPK
0Example and test programs
0      EXAMPLE, JACOBI, TEST, TEST2
1Sec. 2                       MP User's Guide                     Page  8
+                             _______________
02. EXAMPLE PROGRAM
+__________________
0The  following  program  illustrates  the  straightforward use of the MP
 package without the aid of Augment.  An example using Augment  is  given
 in Section 4.3.
0
0C $$                   ******  EXAMPLE  ******
 C
 C THIS PROGRAM COMPUTES PI AND EXP(PI*SQRT(163/9)) TO 100
 C DECIMAL PLACES, AND EXP(PI*SQRT(163)) TO 90 DECIMAL PLACES,
 C AND WRITES THEM ON LOGICAL UNIT 6.
 C
 C CORRECT OUTPUT (EXCLUDING HEADINGS) IS AS FOLLOWS
 C
 C                  3.14159265358979323846264338327950288419716939937510
 C                    58209749445923078164062862089986280348253421170680
 C             640320.00000000060486373504901603947174181881853947577148
 C                    57603665918194652218258286942536340815822646477590
 C 262537412640768743.99999999999925007259719818568887935385633733699086
 C                    2707537410378210647910118607312951181346
 C
 C WORKING SPACE IN BLANK COMMON, PARAMETERS IN
 C COMMON /MPCOM/ (WHICH HAS LENGTH 23)
       COMMON R
       COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY
       INTEGER B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY(16), R(500)
 C
 C WE HAVE T .LE. 62 IF WORDLENGTH AT LEAST 16 BITS, AND WORKING SPACE
 C IS AT MOST 500 WORDS (LESS IF WORDLENGTH IS GREATER THAN 16 BITS).
 C
 C VARIABLES NEED T+2 .LE. 64 WORDS AND ALLOW 110 CHARACTERS FOR
 C DECIMAL OUTPUT
       INTEGER PI(64), X(64), C(110)
 C
 C CALL MPSET2 TO SET OUTPUT LOGICAL UNIT = 6 AND EQUIVALENT
 C NUMBER OF DECIMAL PLACES TO AT LEAST 110.  THE LAST THREE
 C PARAMETERS ARE THE DIMENSIONS OF PI (OR X) AND THE LOWER
 C AND UPPER INDICES OF BLANK COMMON AVAILABLE TO MP.
       CALL MPSET2 (6, 110, 64, 1, 500)
 C
 C COMPUTE MULTIPLE-PRECISION PI
       CALL MPPI(PI)
 C
 C CONVERT TO PRINTABLE FORMAT (F110.100) AND WRITE
       CALL MPOUT (PI, C, 110, 100)
       WRITE (LUN, 10) B, T, C
    10 FORMAT (32H1EXAMPLE OF MP PACKAGE,   BASE =, I9,
      $  12H,   DIGITS =, I4 /// 11H PI TO 100D //
      $  11X, 60A1 / 21X, 50A1)
 C
 C SET X = SQRT(163/9), THEN MULTIPLY BY PI
       CALL MPQPWR (163, 9, 1, 2, X)
       CALL MPMUL (X, PI, X)
1Sec. 2                       MP User's Guide                     Page  9
+                             _______________
0C SET X = EXP(X)
       CALL MPEXP (X, X)
 C
 C CONVERT TO PRINTABLE FORMAT AND WRITE
       CALL MPOUT (X, C, 110, 100)
       WRITE (LUN, 20) C
    20 FORMAT (/ 28H EXP(PI*SQRT(163/9)) TO 100D //
      $        11X, 60A1 / 21X, 50A1)
 C
 C SET X = X**3 = EXP(PI*SQRT(163))
       CALL MPPWR (X, 3, X)
 C
 C WRITE IN FORMAT F110.90
       CALL MPOUT (X, C, 110, 90)
       WRITE (LUN, 30) C
      $        1X, 70A1 / 21X, 40A1)
       WRITE (LUN, 40) MXSPTR
    40 FORMAT (/ 21H END OF EXAMPLE, USED, I4,
      $          23H WORDS OF WORKING SPACE //)
       STOP
       END
1Sec. 3.1                     MP User's Guide                     Page 10
+                             _______________
03. PARAMETERS IN COMMON /MPCOM/
+_______________________________
0
 3.1 Description of the parameters
+_________________________________
0The  labelled  COMMON  block,  COMMON /MPCOM/, contains the following 23
 parameters which are used by various routines  in  the  MP  package  and
 referred  to  in  Section  6.  In the description below, 'default value'
 means the value set by a call to MPSET2 (or MPINIT or MPSET, since  they
 call MPSET2).  All the parameters are of type INTEGER.
0
0(1) BASE or B    The   base    of    the    multiple-precision    number
                  representation.   B  must  be  at least 2, and 8*B**2-1
                  must be representable as  a  single-precision  integer.
                  The   default   value  is  the  largest  power  of  two
                  satisfying this condition.
0(2) NUMDIG or T  The  number  of  multiple-precision  digits (must be at
                  least 2).  The default  value  is  sufficient  to  give
                  accuracy equivalent to at least DECPL decimal places.
0(3) MAXEXP or M  The maximum  exponent  of  multiple-precision  numbers.
                  Exponents  in  the  range 1-M, ... , +M are allowed, so
                  numbers from B**(-M) to almost B**M are  representable.
                  M  must  be  greater  than  T.  The default setting and
                  maximum allowable value is int (MXINT/4).
0(4) LUN          The logical unit for output by routines MPFOUT, MPOUTF,
                  MPDUMP, MPERR, MPERRM etc.  The default  value  is  the
                  first  argument  of  the  last  call to MPSET2 (or 6 if
                  MPINIT as distributed with the MP package is called).
0(5) MXR          Words MNSPTR to MXR of blank COMMON are used as working
                  space  by  the  MP   package.    MXR-MNSPTR   must   be
                  sufficiently large (see Section 6).  The default values
                  of MNSPTR and MXR are the 4-th and  5-th  arguments  of
                  the last call to MPSET2 (or 1 and 500 if MPINIT is used
                  as distributed with the MP package).
0(6) SPTR         The  current  pointer to the top of the working storage
                  stack  (plus  one).    SPTR  should  be  in  the  range
                  MNSPTR <= SPTR <= MXR + 1.  The default initial setting
                  is MNSPTR.
0(7) MXSPTR       The maximum value of SPTR during the current run.   The
                  default initial setting is MNSPTR.
0(8) MNSPTR       The minimum value of SPTR during the current run.   See
                  the description of MXR above.
0(9) MXEXPN       The maximum exponent of any  multiple-precision  number
                  during the current run.  The default initial setting is
                  -M.  Set to M+1 if multiple-precision overflow  occurs.
1Sec. 3.1                     MP User's Guide                     Page 11
+                             _______________
0(10) MNEXPN      The  minimum  exponent of any multiple-precision number
                  during the current run.  The default initial setting is
                  M+1.   Set  to  -M  if  a  multiple-precision underflow
                  occurs.
0(11) RNDRL       Integer  in the range 0, ... , 3 indicating the type of
                  rounding to be performed.  Generally,
0                   0 means truncated (chopped) arithmetic,
                    1 means rounded arithmetic,
                    2 means round down (towards -infinity),
                    3 means round up (towards +infinity).
0                 However,  these  interpretations  do  not always apply.
                  For details of the effect of RNDRL, see Section 6.  The
                  default setting is 0 (truncated arithmetic).
0(12) KTUNFL      Counter for the number of multiple-precision underflows
                  which  have  occurred.   The default initial setting is
                  zero.
0(13) MXUNFL      If  positive,  execution  is  terminated  after  MXUNFL
                  multiple-precision underflows have occurred.  If  zero,
                  any number of multiple-precision underflows are allowed
                  (this is the default).
0(14) DECPL       A  conservative  estimate of the 'equivalent' number of
                  floating-point decimal places.  The default setting  is
                  the  second  argument in the last call to MPSET2 (or 43
                  if  MPINIT  is  used  as  distributed).   The   default
                  settings of B and T satisfy
0                   B**(T - 1) >= 10**(DECPL - 1)
0(15) MT2         The  dimension  of  arrays  used for multiple-precision
                  variables.  The default setting is the  third  argument
                  in  the last call to MPSET2 (or 27 if MPINIT is used as
                  distributed).
0(16) MXINT       The  largest  single-precision   integer  of  the  form
                  2**k - 1 such that INTEGER arithmetic with operands and
                  results in the range -MXINT, ... , +MXINT is  performed
                  exactly.   Usually  k  = (wordlength in bits) - 1.  See
                  also the description of MPLARG in Section 6.
0(17) EXWID       The  number  of  characters  in  the exponent field for
                  output of multiple-precision numbers using MPFOUT.  The
                  default value is 6. (Since one position is taken by the
                  exponent character EXPCH and  one  by  the  sign,  this
                  means  that  exponents  in the range -9999 to +9999 are
                  allowable with the default setting.)
0(18) INRECL      The  input  record  length  assumed by MPFIN.  Must not
                  exceed 80 (the default value) unless  MPFIN  and  MPCHK
                  are modified.
1Sec. 3.2                     MP User's Guide                     Page 12
+                             _______________
0(19) INBASE      The base of 'decimal' numbers read by MPIN, MPFIN  etc.
                  (e.g.  INBASE  = 8 means octal input, INBASE = 16 means
                  hexadecimal input).  Must be in the range 2, ... ,  16.
                  The default setting is ten, for decimal input.
0(20) OUTBAS      The base for 'decimal' output  by  MPOUT,  MPFOUT  etc.
                  Default  setting ten, must be in the range 2, ... , 16.
0(21) EXPCH       The exponent character used by MPFOUT.  Must not  be  a
                  digit, sign or blank.  The default setting is 'E'.
0(22) CHWORD      The number of characters per  single-precision  integer
                  word.  The default setting is determined by MPUPW.
0(23) ONESCP      1 or 0 if the machine  uses  one's  complement  integer
                  arithmetic,   0  otherwise.   The  default  setting  is
                  determined by MPSET2 and MPUPW so that MPUPW works.
03.2 Methods of setting the parameters
+_____________________________________
0It would, of course, be possible to declare the parameters B, T, M, LUN,
 ...  ,  ONESCP  in  COMMON  /MPCOM/  and  initialize  them by assignment
 statements or by a BLOCK DATA subprogram.  However, this method  is  not
 recommended.   Instead, we recommend calling MPINIT or MPSET2 to set the
 default values (see above and Section 6), and then using MPPARB or PARAM
 to change the setting of any parameters for which the default values are
 not desired.
0For  example,  to  use  the  MP  package  with 20 place rounded floating
 decimal arithmetic, we could do the following -
0      Without using Augment               Using Augment
+      _____________________               _____________
0      CALL MPINIT (MP)                    INITIALIZE MP
       CALL MPPARB (20, 'NUMDIG')          PARAM ('NUMDIG') = 20
       CALL MPPARB (10, 'BASE')            PARAM ('BASE') = 10
       CALL MPPARB (1, 'RNDRL')            PARAM ('RNDRL') = 1
0MPPARC may be used instead of MPPARB.  It is less mnemonic, but slightly
 faster, so might be preferred inside a loop. For details see Section  6.
0The  current  setting  of  the parameters may be determined by MPPARA or
 PARAM.  For example, MPPARA ('NUMDIG') will return the current number of
 multiple-precision  digits.   For  users of Augment, PARAM ('NUMDIG') on
 the right-hand side of an assignment  statement  will  return  the  same
 information.
0Augment users may prefer to set the first  three  parameters  in  COMMON
 /MPCOM/  with  the  field  functions  BASE, NUMDIG and MAXEXP instead of
 PARAM.  For example, if X is a variable of type MULTIPLE, the number  of
 multiple-precision digits may be set to 20 by the statement
0      NUMDIG (X) = 20
0For further details, see Section 5.2.
1Sec. 3.3                     MP User's Guide                     Page 13
+                             _______________
03.3 Efficiency considerations
+_____________________________
0For efficiency choose B fairly large, subject to the restrictions  given
 in Section 1.2.  For example, if the wordlength is
0      48 bits, could use B = 4194304 or 1000000,
       36 bits, could use B = 65536 or 100000,
       32 bits, could use B = 16384 or 10000,
       24 bits, could use B = 1024 or 1000,
       18 bits, could use B = 128 or 100,
       16 bits, could use B = 64 or 10,
       12 bits, could use B = 16 or 10.
0Of  course, for special purposes such as simulation of decimal or binary
 arithmetic, B may be set as small as 2.   If  a  significant  amount  of
 decimal  input/output  is  to be performed, it may be best to choose B a
 power of ten (MPIN and MPOUT will take advantage of this).
0Avoid  multiplication  and  division  by  multiple-precision numbers, as
 these take O(T**2) operations, whereas multiplication  and  division  by
 (single-precision)  integers take O(T) operations.  See the descriptions
 of MPDIV, MPDIVI, MPMUL and MPMULI in Section 6.
0The  default  setting of RNDRL is 0 (i.e. truncated arithmetic). The use
 of other values of RNDRL is more expensive in both time and space.
1Sec. 4.1                     MP User's Guide                     Page 14
+                             _______________
04. THE AUGMENT INTERFACE
+________________________
04.1 Introduction
+________________
0Augment  is  a preprocessor which allows the introduction of nonstandard
 types (e.g. multiple-precision variables) into  Fortran  programs.   For
 details,  see  -  F.D.  Crary,  'A versatile precompiler for nonstandard
 arithmetics', ACM Trans. Math. Software 5 (June 1979), 204-217,  and the
 references given there.
0A 'description deck' has been written to enable Augment to  be  used  in
 conjunction  with  the  MP package.  This greatly simplifies the task of
 writing a program for a multiple-precision computation, or converting  a
 single (or double) precision routine to multiple precision.
0For example, if Augment is used  we  can  declare  X,  Y  and  Z  to  be
 multiple-precision variables by a declaration such as
0      MULTIPLE X, Y, Z
 or
       IMPLICIT MULTIPLE (A-H, O-Z)
0and subsequently write expressions such as
0       X = Y + Z*EXP(X+1)/Y
0instead of the equivalent
0       CALL MPADDI (X, 1, TEMP)
        CALL MPEXP (TEMP, TEMP)
        CALL MPMUL (Z, TEMP, TEMP)
        CALL MPDIV (TEMP, Y, TEMP)
        CALL MPADD (Y, TEMP, X)
0The Augment interface can be used  with  MP  version  780420,  or  later
 versions.   For  more details, see - R. P. Brent, J. A. Hooper and J. M.
 Yohe, 'An Augment interface for  Brent's  multiple-precision  arithmetic
 package', ACM Trans. Math. Software 6 (1980), 146-149.
0A convenient summary of the operations implemented in the MP package and
 the methods of invoking them via Augment is given in Section 5.2.
1Sec. 4.2                     MP User's Guide                     Page 15
+                             _______________
04.2 The description deck
+________________________
0The 'description deck' which describes the MP package is as follows.
0COMMENT  AUGMENT DESCRIPTION DECK FOR THE MULTIPLE-PRECISION ARITHMETIC
          PACKAGE OF R. P. BRENT.
0         THREE TYPES OF VARIABLE ARE DEFINED HERE -
0               MULTIPLE   (STANDARD MULTIPLE-PRECISION NUMBERS),
                MULTIPAK   (PACKED MULTIPLE-PRECISION NUMBERS), AND
                INITIALIZE (USED ONLY AS A DEVICE TO PERSUADE
                           AUGMENT TO INITIALIZE THE MP PACKAGE).
0         WORKING SPACE SHOULD BE ALLOCATED AND THE MP PACKAGE
          INITIALIZED BY THE DECLARATION
0               INITIALIZE MP
0         IN THE MAIN PROGRAM.
0         THIS DESCRIPTION DECK ASSUMES THAT EACH MP NUMBER REQUIRES NO
          MORE THAN 27 WORDS (14 IN PACKED FORMAT).  THIS IS SUFFICIENT
          FOR ABOUT 43 DECIMAL PLACE ACCURACY ON A 16-BIT MACHINE, AND
          HIGHER ACCURACY ON A MACHINE WITH A LONGER WORDLENGTH.
0         SEE COMMENTS IN ROUTINE MPINIT FOR THE METHOD OF CHANGING THE
          PRECISION, ALSO REGARDING COMMON DECLARATIONS.
0*DESCRIBE MULTIPAK
0DECLARE INTEGER (14),  KIND SAFE SUBROUTINE, PREFIX MPK
0CONVERSION CTP (CIM, INTEGER, $, UPWARD),  CTP (CAM, HOLLERITH)
0SERVICE COPY (=MPSTR)
0*DESCRIBE MULTIPLE
0DECLARE INTEGER (27),  KIND SAFE SUBROUTINE, PREFIX MP
0OPERATOR + (, NULL UNARY, PRV, $),  - (NEG, UNARY),
          + (ADD, BINARY3, PRV, $, $, $, COMM),  * (MUL),
          - (SUB,,,,,, NONCOMM),  / (DIV),  ** (PWR2),
0         + (ADDI,,,, INTEGER),  * (MULI),  / (DIVI),  ** (PWR),
          * (IMUL,,, INTEGER, $, $, NONCOMM),
0         .EQ. (EQ, BINARY2, PRV, $, LOGICAL, COMM),  .NE. (NE),
          .GE. (GE,,,,, NONCOMM),  .GT. (GT),  .LE. (LE),  .LT. (LT),
0         + (, NULL UNARY, PRV, MULTIPAK),
          - (NEG, UNARY, PRV, MULTIPAK),
          + (KADD, BINARY3, PRV, MULTIPAK, MULTIPAK, $, COMM),
          * (KMUL),  - (KSUB,,,,,, NONCOMM),  / (KDIV),  ** (PWR2),
1Sec. 4.2                     MP User's Guide                     Page 16
+                             _______________
0         * (KMLI,,,, INTEGER),  / (KDVI),  ** (PWR),
          * (KIML,,, INTEGER, MULTIPAK, $, NONCOMM),
0         .EQ. (EQ, BINARY2, PRV, MULTIPAK, LOGICAL, COMM),  .NE. (NE),
          .GE. (GE,,,,, NONCOMM),  .GT. (GT),  .LE. (LE),  .LT. (LT)
0TEST     MPSIGA (SIGA, INTEGER)
0
 FIELD    SGN (SIGA, SIGB,  ($),  INTEGER),  EXPON (EXPA, EXPB),
          BASE (BASA, BASB),  NUMDIG (DIGA, DIGB),  MAXEXP (MEXA, MEXB),
          DIGIT (DGA, DGB,  ($, INTEGER)),
          PARAM (PARA, PARB,  (HOLLERITH)),  FIO (FIN, FOUT,  (INTEGER),
          $),  UNFIO (UNFR, UNFW),
          EXPON (EXPA, EXPB, (MULTIPAK),  INTEGER),  NUMDIG (DIGA, DIGB),
          SGN (SIGA, SIGB)
0
 FUNCTION ABS (ABS,  ($),  $),  DABS (ABS),  ASIN (ASIN),  DASIN (ASIN),
          ARSIN (ASIN),  DARSIN (ASIN),  ARCSIN (ASIN),  ATAN (ATAN),
          DATAN (ATAN),  ARTAN (ATAN),  DARTAN (ATAN),  ARCTAN (ATAN),
          CMF (CMF),  CEIL (CEIL),  FLOOR (FLOR),  CMIM (CMIM),
          COS (COS),  DCOS (COS),  COSH (COSH),  DCOSH (COSH),
          DAW (DAW),  DDAW (DAW),  EI (EI),  DEI (EI),  ERF (ERF),
          DERF (ERF),  ERFC (ERFC),  DERFC (ERFC),  EXP (EXP),
          DEXP (EXP),  EXP1 (EXP1),  FRAC (CMF),  GAM (GAM),
          DGAM (GAM),  GAMMA (GAM),  DGAMMA (GAM),  AINT (CMIM),
          DINT (CMIM),  LI (LI),  DLI (LI),  LN (LN),  LOG (LN),
          ALOG (LN),  DLOG (LN),  LOG10 (LG10),  ALOG10 (LG10),
          DLOG10 (LG10),  LNGM (LNGM),  ALNGM (LNGM),  DLNGM (LNGM),
          LNGS (LNGS),  LNS (LNS),  REC (REC),  SIN (SIN),  DSIN (SIN),
          SINH (SINH),  DSINH (SINH),  SQRT (SQRT),  DSQRT (SQRT),
          TAN (TAN),  DTAN (TAN),  TANH (TANH),  DTANH (TANH),
0         SNGL (STR),  DBLE (STR),
0         ART1 (ART1,  (INTEGER)),  LN (LNI),  LNI (LNI),  LOG (LNI),
          ALOG (LNI),  DLOG (LNI),  ZETA (ZETA),
0         CAM (CAM),  CAM (CAM,  (HOLLERITH)),
0         MAX (MAX,  ($, $)),  AMAX1 (MAX),  DMAX1 (MAX),  MIN (MIN),
          AMIN1 (MIN),  DMIN1 (MIN),  DIM (DIM),  DDIM (DIM),
          GCD (GCDA),  ATAN2 (ATN2),  DATAN2 (ATN2),  ARTAN2 (ATN2),
          MOD (MOD),  AMOD (MOD),  DMOD (MOD),
          SIGN (SIGN),  DSIGN (SIGN),
0         BESJ (BESJ,  ($, INTEGER)),  ROOT (ROOT),
0         MPINF (INF (SUBROUTINE),   ($, INTEGER, INTEGER, HOLLERITH),
          LOGICAL),  MPOUTF (OUTF (SUBROUTINE)),
          MPINF (INF (SUBROUTINE),   ($, INTEGER, INTEGER, INTEGER)),
          MPOUTF (OUTF (SUBROUTINE)),
0         INTEGR (INTG,  ($),  LOGICAL),
1Sec. 4.2                     MP User's Guide                     Page 17
+                             _______________
0         COMP (COMP,  ($, $),  INTEGER),  CMPA (CMPA),
          COMP (CMPI,  ($, INTEGER)),  COMP (CMPR,  ($, REAL)),
          COMP (CMPD,  ($, DOUBLE PRECISION)),
          COMP (CMPQ,  ($, INTEGER, INTEGER)),
0         ADDQ (ADDQ,  ($, INTEGER, INTEGER),  $),  MULQ (MULQ),
          QPWR (QPWR,  (INTEGER, INTEGER, INTEGER, INTEGER)),
          CQM (CQM,  (INTEGER, INTEGER)),  CTM (CQM),
0         GAM (GAMQ),  DGAM (GAMQ),  GAMQ (GAMQ),
0         BERN (BERN,  (INTEGER, INTEGER),  MULTIPAK),
0         INT (CMI (SUBROUTINE),   ($),  INTEGER),
          IDINT (CMI (SUBROUTINE)),  IFIX (CMI (SUBROUTINE)),
0         ABS (ABS, (MULTIPAK), MULTIPAK), ATAN (ATAN,, $), COS (COS),
          COSH (COSH),  EXP (EXP),  LN (LN),  LOG (LN),  ALOG (LN),
          LOG10 (LG10),  ALOG10 (LG10),  SIN (SIN),  SINH (SINH),
          SQRT (SQRT),  TAN (TAN),  TANH (TANH),
          EXP1 (EXP1), LNS (LNS),
0         MAX (MAX,  (MULTIPAK, MULTIPAK)),  MIN (MIN),  DIM (DIM),
          ATAN2 (ATN2),  MOD (MOD),  SIGN (SIGN),
          ROOT (ROOT, (MULTIPAK, INTEGER))
0CONVERSION CTM (CDM, DOUBLE PRECISION, $, UPWARD),  CTM (CIM, INTEGER),
          CTM (CRM, REAL),  CTM (UNPK, MULTIPAK),  CTM (CAM, HOLLERITH),
          CTD (CMD (SUBROUTINE),  $, DOUBLE PRECISION, DOWNWARD),
          CTI (CMI (SUBROUTINE),, INTEGER),  CTR (CMR (SUBROUTINE),,
          REAL),  CTP (PACK,, MULTIPAK)
0SERVICE COPY (STR)
0*DESCRIBE INITIALIZE
0DECLARE INTEGER (1),  KIND SAFE SUBROUTINE, PREFIX MPI
 SERVICE COPY (STR),  INITIAL (NIT)
0COMMENT    END OF AUGMENT DESCRIPTION DECK FOR MP PACKAGE
1Sec. 4.3                     MP User's Guide                     Page 18
+                             _______________
04.3 JACOBI - an example of the use of Augment
+_____________________________________________
0The program which follows illustrates the use of the MP package with the
 Augment  interface.   It  is  quite straightforward, and does not by any
 means  illustrate  all  the  possibilities  provided  by   the   Augment
 interface.
0
 < Machine-dependent statement(s) to execute Augment >
 PRINT SUPPRESS
 < Machine-dependent statement(s) to insert the description deck >
0*BEGIN
 *DISABLE PRINT, OUTPUT
 *ENABLE SOURCE
0< Machine-dependent statement(s) to (later) execute Fortran compiler >
0C
 C PROGRAM TO VERIFY AN IDENTITY OF JACOBI USING THE MP
 C PACKAGE AND AUGMENT.
 C
 C THE PROGRAM READS A NUMBER X IN FREE-FIELD FORMAT ACCEPTABLE TO
 C MPIN.  IF X IS NON-POSITIVE IT HALTS.  OTHERWISE IT COMPUTES
 C AND PRINTS FN(X), FN(1/X) AND (FN(X)-FN(1/X))/FN(X),
 C WHERE  FN(X) IS THE SUM FROM N = -INFINITY TO +INFINITY OF
 C SQRT(X)*EXP(-PI*(N*X)**2).
 C THE IDENTITY VERIFIED IS
 C                           FN(X) = FN(1/X)
 C
 C DECLARE VARIABLES AND INITIALIZE MP PACKAGE.  ON SOME SYSTEMS BLANK
 C COMMON MUST BE DECLARED HERE - SEE COMMENTS IN ROUTINE MPINIT.
 C
       MULTIPLE X, F1, F2, FN
       INITIALIZE MP
 C
 C SET INPUT RECORD LENGTH TO 72 (TO AVOID READING
 C SEQUENCE NUMBERS), DEFAULT VALUE IS 80.
 C
       PARAM (6HINRECL) = 72
 C
 C READ MP X FROM UNIT 5 IN FREE FORMAT, STOP IF X NOT POSITIVE.
 C
    10 X = FIO (5)
       IF (X.LE.0) STOP
 C
 C WRITE HEADING, X, FN(X), AND FN(1/X) TO 40S IN STANDARD FORMAT
 C
       WRITE (6, 20)
    20 FORMAT (//41H X, FN(X), FN(1/X), (FN(X)-FN(1/X))/FN(X)/)
       FIO (40) = X
       F1 = FN(X)
       FIO (40) = F1
       F2 = FN(1/X)
       FIO (40) = F2
1Sec. 4.3                     MP User's Guide                     Page 19
+                             _______________
0C WRITE (F1-F2)/F1 TO 6 SIGNIFICANT FIGURES.
 C NOTE THAT AN MP EXPRESSION CAN BE WRITTEN.
 C
       FIO (6) = (F1 - F2)/F1
       GO TO 10
       END
 C
 < Machine-dependent statement(s) to (later) execute Fortran compiler >
 C
       FUNCTION FN(X)
 C
 C RETURNS FN(X) = THE SUM FROM N = -INFINITY TO +INFINITY OF
 C SQRT(X)*EXP(-PI*(N*X)**2), ASSUMING X POSITIVE.
 C USES THE OBVIOUS METHOD, SO SLOW IF X SMALL.
 C NOTE THAT X AND FN ARE BOTH TYPE MULTIPLE.
 C TO ILLUSTRATE MIXED-MODE ARITHMETIC, SOME LOCAL
 C VARIABLES ARE DECLARED TO BE PACKED.
 C
       MULTIPLE FN, X
       MULTIPAK TM, FAC, PR
       IF (X.LE.0) CALL MPERR
       FN = 0
 C
 C AUGMENT CAN DEAL WITH THE FOLLOWING EXPRESSION AS IT KNOWS THAT X
 C IS TYPE MULTIPLE, SO CALLS MPCAM TO CONVERT 'PI' TO MULTIPLE.
 C
       TM = EXP(-2HPI*X*X)
       PR = TM
       FAC = TM**2
 C
 C LOOP TO SUM INFINITE SERIES
 C WARNING - NUMBER OF ITERATIONS IS PROPORTIONAL TO 1/X
 C
    10 FN = FN + TM
       PR = PR*FAC
       TM = TM*PR
 C
 C TEST FOR CONVERGENCE BY COMPARING EXPONENTS OF FN AND TM.
 C WE COULD ALSO HAVE SAVED THE OLD VALUE OF FN AND SEEN IF
 C STATEMENT 10 CHANGED IT.
 C
       IF (EXPON(FN)-EXPON(TM).LT.NUMDIG(X)) GO TO 10
       FN = SQRT(X)*(2*FN+1)
       RETURN
       END
 *END
0<Machine-dependent statements to add Augment output file to
  the runstream, link compiled Fortran routines with precompiled
  library of MP routines, and execute with the following data >
 .5
 .3
 1.+1
 1.234567890123456789012345678901234567890123456789
 0
1Sec. 5.1                     MP User's Guide                     Page 20
+                             _______________
05. SUMMARY FOR AUGMENT USERS
+____________________________
05.1 Notes
+_________
0We summarize as briefly as possible how MP routines may be invoked using
 the Augment interface.  Under 'Invocation' we give the  typical  way  of
 invoking  the  operation  described  under  'Description'.   For further
 details of calling sequences etc., the column 'Ref' gives a reference to
 the relevant MP routine, and descriptions of these routines may be found
 in Section 6.  The column headed 'Res' gives the type of the  result  of
 each operation.
0Types are indicated by one-letter abbreviations,
0      D = double-precision real (DOUBLE PRECISION),
       G = packed or unpacked multiple (MULTIPAK or MULTIPLE)
           (indicated by 'MP*' in the column headed 'Description'),
       H = packed Hollerith (as in 'ABCDEF'),
       I = single-precision integer (INTEGER),
       L = Boolean (LOGICAL),
       M = unpacked multiple (MULTIPLE),
       P = packed multiple (MULTIPAK),
       R = single-precision real (REAL).
0Arguments are written as XA, XB, XC, ... , where X indicates  the  type,
 as  above.   XDUMMY  means  a dummy argument of the type indicated by X.
 Such arguments serve no useful purpose, but Augment expects them.
0The   conversion  routines  may  be  invoked  implicitly  by  assignment
 statements, for example -
0      MR = IA
0will be translated by Augment into a call to MPCIM.
0Packed   (i.e.   type   MULTIPAK)  and  unpacked  (i.e.  type  MULTIPLE)
 multiple-precision  variables  may  be  mixed   in   simple   arithmetic
 expressions. For example, in the statement -
0      A = B + C/(D - E*F)
0any  combination  of packed and unpacked multiple-precision variables is
 permissible.  Augment will generate calls on the appropriate  conversion
 routines.  We do not recommend using packed multiple-precision variables
 in  expressions  involving  function  calls,   without   explicit   type
 conversion  (using  CTM  and/or  CTP),  or  unless  it  is  mentioned as
 permissible in the summary below.   Augment  will  not  always  generate
 calls on the appropriate conversion routines.
0Field functions may occur on either side of an assignment statement. For
 example -
0      PARAM ('MXUNFL') = 10
0would set the maximum allowable number of MP underflows to 10.
1Sec. 5.2                     MP User's Guide                     Page 21
+                             _______________
05.2 Summary
+___________
0             Description                    Res Invocation         Ref
+             ___________                    ___ __________         ___
0  Addition
+  ________
0      Sum of two MP numbers                  M  MA + MB           MPADD
       Sum of two packed MP numbers           M  PA + PB           MPKADD
       Sum of MP number and an integer        M  MA + IB           MPADDI
       Sum of MP number and rational IB/IC    M  ADDQ (MA, IB, IC) MPADDQ
0  Division
+  ________
0      Quotient of two MP numbers             M  MA / MB           MPDIV
       Quotient of two packed MP numbers      M  PA / PB           MPKDIV
       Quotient of MP number and an integer   M  MA / IB           MPDIVI
       Quotient of packed MP and integer      M  PA / IB           MPKDVI
0  Multiplication
+  ______________
0      Product of two MP numbers              M  MA * MB           MPMUL
       Product of two packed MP numbers       M  PA * PB           MPKMUL
       Product of MP number and integer       M  MA * IB           MPMULI
       Product of packed MP and integer       M  PA * IB           MPKMLI
       Product of integer and MP number       M  IA * MB           MPIMUL
       Product of integer and packed MP       M  IA * PB           MPKIML
       Product of MP and rational IB/IC       M  MULQ (MA, IB, IC) MPMULQ
0  Reciprocal
+  __________
0      Reciprocal of MP number                M  REC (MA)          MPREC
0  Subtraction
+  ___________
0      Difference of two MP numbers           M  MA - MB           MPSUB
       Difference of two packed MP numbers    M  PA - PB           MPKSUB
0  Powers and roots
+  ________________
0      Raise MP* number to integer power      M  GA ** IB          MPPWR
       Raise MP number to MP power            M  MA ** MB          MPPWR2
       Raise packed MP to packed MP power     M  PA ** PB          MPPWR2
       Raise rat. IA/IB to rat. power IC/ID   M  QPWR(IA,IB,IC,ID) MPQPWR
       IB-th root of MP* number               M  ROOT (GA, IB)     MPROOT
0
   Elementary functions
+  ____________________
0      Arcsin of MP number (in radians)       M  ASIN (MA)         MPASIN
       Arctan of MP* number (in radians)      M  ATAN (GA)         MPATAN
       Arctan of 1/IA, IA > 1                 M  ART1 (IA)         MPART1
       Arctan of MA/MB (both MP numbers)      M  ATAN2 (MA, MB)    MPATN2
       Arctan of PA/PB (both packed MP)       M  ATAN2 (PA, PB)    MPATN2
       Cosine of MP* number (radians)         M  COS (GA)          MPCOS
       Hyperbolic cosine of MP* number        M  COSH (GA)         MPCOSH
1Sec. 5.2                     MP User's Guide                     Page 22
+                             _______________
0             Description                    Res Invocation         Ref
+             ___________                    ___ __________         ___
0  Elementary functions continued
+  ______________________________
0      Exponential of MP* number              M  EXP (GA)          MPEXP
       (Exponential - 1) of MP* number        M  EXP1 (GA)         MPEXP1
       Natural logarithm of MP* number        M  LOG (GA)          MPLN
       Logarithm (base 10) of MP* number      M  LOG10 (GA)        MPLG10
       Natural logarithm of (1 + MP* number)  M  LNS (GA)          MPLNS
       Natural logarithm of integer > 0       M  LOG (IA)          MPLNI
       Sine of MP* number (radians)           M  SIN (GA)          MPSIN
       Hyperbolic sine of MP* number          M  SINH (GA)         MPSINH
       Square root of MP* number              M  SQRT (GA)         MPSQRT
       Tangent of MP* number (radians)        M  TAN (GA)          MPTAN
       Hyperbolic tangent of MP* number       M  TANH (GA)         MPTANH
0
   Special functions
+  _________________
0      Bessel function of first kind,
          MP argument MA, integer order IB    M  BESJ (MA, IB)     MPBESJ
       Dawson's integral of MP argument       M  DAW (MA)          MPDAW
       Exponential integral of MP argument    M  EI (MA)           MPEI
       Error function of MP number            M  ERF (MA)          MPERF
       Complementary error function           M  ERFC (MA)         MPERFC
       Gamma function of MP number            M  GAM (MA)          MPGAM
       Gamma function of rational IA/IB       M  GAM (IA, IB)      MPGAMQ
       GCD of MP numbers which are integers   M  GCD (MA, MB)      MPGCDA
       Logarithmic integral of MP number      M  LI (MA)           MPLI
       Logarithm of Gamma function            M  LNGM (MA)         MPLNGM
       Riemann zeta function of integer > 1   M  ZETA (IA)         MPZETA
0
   Constants
+  _________
0      Bernoulli numbers (result is array)    P  BERN (IA, IB)     MPBERN
       MP machine precision                   M  CTM ('EPS')       MPEPS
       Euler's constant                       M  CTM ('EUL')       MPEUL
       Largest positive MP number             M  CTM ('MAXR')      MPMAXR
       Smallest positive MP number            M  CTM ('MINR')      MPMINR
       Pi                                     M  CTM ('PI')        MPPI
       Decimal constants, e.g. 6.3E5          M  CTM ('6.3E5$')    MPCAM
0
   Input and output
+  ________________
0      Read MP number from unit IA in
          free format                         M  FIO (IA)          MPFIN
0      Write MP number MB with IA decimal
          places on default output unit          FIO (IA) = MB     MPFOUT
0      Read IB words from unit IC under
          format HD and convert to MP         L  MPINF(MA,IB,IC,HD)
          number MA, result error flag                             MPINF
1Sec. 5.2                     MP User's Guide                     Page 23
+                             _______________
0             Description                    Res Invocation         Ref
+             ___________                    ___ __________         ___
0  Input and output continued
+  __________________________
0      Write IB words representing MP
          number MA on default unit, with
          IC places after decimal point,      L  MPOUTF(MA,IB,IC,HD)
          under format HD, result error flag                       MPOUTF
0      Read MP number from unit IA,
          unformatted                         M  UNFIO (IA)        MPUNFR
0      Write MP number MB on unit IA,
          unformatted                            UNFIO (IA) = MB   MPUNFW
0
0  Conversion
+  __________
0      Double-precision to multiple           M  CTM (DA)          MPCDM
       Integer to multiple                    M  CTM (IA)          MPCIM
       Real to multiple                       M  CTM (RA)          MPCRM
       Rational IA/IB to multiple             M  CTM (IA, IB)      MPCQM
       Packed multiple to unpacked multiple   M  CTM (PA)          MPUNPK
       Packed Hollerith to multiple           M  CTM (HA)          MPCAM
       Packed Hollerith stored in integer
          or integer array to multiple        M  CAM (IA)          MPCAM
       Multiple to double-precision           D  CTD (MA)          MPCMD
       Multiple to integer (truncated)        I  CTI (MA)          MPCMI
       Multiple to real                       R  CTR (MA)          MPCMR
       Unpacked multiple to packed multiple   P  CTP (MA)          MPPACK
0
0  Comparison (result is +1, 0 or -1)
+  __________
0      Compare absolute value of MP numbers   I  CMPA (MA, MB)     MPCMPA
       Compare MP number with integer         I  COMP (MA, IB)     MPCMPI
       Compare MP number with real            I  COMP (MA, RB)     MPCMPR
       Compare MP number with double-prec.    I  COMP (MA, DB)     MPCMPD
       Compare MP number with rat. IB/IC      I  COMP (MA, IB, IC) MPCMPQ
       Compare two MP numbers                 I  COMP (MA, MB)     MPCOMP
0
0  Relational
+  __________
0      Compare MP* numbers for equality       L  GA .EQ. GB        MPEQ
          (Similarly for .GE.,                L  GA .GE. GB        MPGE
          .GT., .LE., .LT. and .NE.)                etc.           etc.
0
   Test
+  ____
0      Three-way branch (n1, n2, n3 labels)      IF (MA) n1,n2,n3  MPSIGA
1Sec. 5.2                     MP User's Guide                     Page 24
+                             _______________
0             Description                    Res Invocation         Ref
+             ___________                    ___ __________         ___
0  Field functions
+  _______________
0      Sign of MP* number                     I  SGN (GA)          MPSIGA
                                                                   MPSIGB
       Exponent of MP* number                 I  EXPON (GA)        MPEXPA
                                                                   MPEXPB
       IB-th digit of MP number               I  DIGIT (MA, IB)    MPDGA
                                                                   MPDGB
       Number of MP digits (T)                I  NUMDIG (GDUMMY)   MPDIGA
                                                                   MPDIGB
       Maximum exponent of MP numbers (M)     I  MAXEXP (MDUMMY)   MPMEXA
                                                                   MPMEXB
       Multiple-precision base (B)            I  BASE (MDUMMY)     MPBASA
                                                                   MPBASB
       Parameter in COMMON /MPCOM/,           I  PARAM (HA)        MPPARA
          argument HA may be any of 'BASE',                        MPPARB
          'NUMDIG', 'MAXEXP', 'LUN', 'SPTR',
          'MXSPTR', 'MNSPTR', 'MXEXPN',
          'MNEXPN', 'RNDRL', 'KTUNFL',
          'MXUNFL', 'DECPL', 'MT2', 'MXINT',
          'EXWID', 'INRECL', 'INBASE',
          'OUTBAS', 'EXPCH', 'CHWORD',
          or 'ONESCP'
0      See under 'Input and output'           M  FIO (IA)          MPFIN
                                                                   MPFOUT
       See under 'Input and output'           M  UNFIO (IA)        MPUNFR
                                                                   MPUNFW
0  Miscellaneous
+  _____________
0      Unary minus of MP number               M  -MA               MPNEG
       Unary minus of packed MP number        P  -PA               MPNEG
       Assignment of MP number (result = MA)  M  MA                MPSTR
       Assignment of packed MP number         P  PA                MPSTR
       Absolute value of MP number            M  ABS (MA)          MPABS
       Absolute value of packed MP number     P  ABS (PA)          MPABS
       Fractional part of MP number           M  FRAC (MA)         MPCMF
       Integer part of MP number              M  INT (MA)          MPCMIM
       Ceiling function of MP number          M  CEIL (MA)         MPCEIL
       Floor function of MP number            M  FLOOR (MA)        MPFLOR
       Returns true if MA is exact integer    L  INTEGR (MA)       MPINTG
       Maximum of two MP numbers              M  MAX (MA, MB)      MPMAX
       Maximum of two packed MP numbers       M  MAX (PA, PB)      MPMAX
       Minimum of two MP numbers              M  MIN (MA, MB)      MPMIN
       Minimum of two packed MP numbers       M  MIN (PA, PB)      MPMIN
       Max (0, MA - MB) for MP numbers        M  DIM (MA, MB)      MPDIM
       Max (0, PA - PB) for packed numbers    M  DIM (PA, PB)      MPDIM
       MA - int(MA/MB)*MB for MP numbers      M  MOD (MA, MB)      MPMOD
       PA - int(PA/PB)*PB for packed numbers  M  MOD (PA, PB)      MPMOD
       abs(MA)*sign(MB) for MP numbers        M  SIGN (MA, MB)     MPSIGN
       abs(PA)*sign(PB) for packed numbers    M  SIGN (PA, PB)     MPSIGN
       Initialization of MP package              INITIALIZE MP     MPINIT
1Sec. 5.3                     MP User's Guide                     Page 25
+                             _______________
05.3 Synonyms
+____________
0In some cases Augment will recognize synonyms for the names given in the
 column headed 'Invocation' in the summary above.  For example, AMAX1 may
 be  used  in  place  of MAX.  A list of synonyms follows.  In most cases
 these are usable only for unpacked multiple-precision arguments.  It  is
 easy  to  modify or expand the list by making appropriate changes to the
 Augment description deck.
0      Generic name          Synonyms
+      ____________          ________
0         ABS                DABS
          AINT               DINT
          ASIN               ARCSIN, ARSIN, DARSIN, DASIN
          ATAN               ARCTAN, ARTAN, DARTAN, DATAN
          ATAN2              ARTAN2, DATAN2
          COS                DCOS
          COSH               DCOSH
          DAW                DDAW
          DIM                DDIM
          EI                 DEI
          ERF                DERF
          ERFC               DERFC
          EXP                DEXP
          GAM                DGAM, DGAMMA, GAMMA,
                             GAMQ (rational argument only)
          INT                IDINT, IFIX
          LI                 DLI
          LOG                ALOG, DLOG, LN, LNI (integer argument only)
          LNGM               ALNGM, DLNGM
          LOG10              ALOG10, DLOG10
          MAX                AMAX1, DMAX1
          MIN                AMIN1, DMIN1
          SIGN               DSIGN
          MOD                AMOD, DMOD
          SIN                DSIN
          SINH               DSINH
          SQRT               DSQRT
          TAN                DTAN
          TANH               DTANH
05.4 Reserved words
+__________________
0When writing programs which use MP via the Augment interface, it is best
 to  avoid  using  the  following  identifiers except with their reserved
 meanings as indicated above -
0      INITIALIZE,
       MPxxxx (for any letters or digits xxxx),
       MULTIPAK,
       MULTIPLE.
0Care  should  also be taken to avoid using the names given in the column
 headed 'Invocation' in the summary above, or in  the  list  of  synonyms
 above, in a manner which would confuse Augment.
1Sec. 6                       MP User's Guide                     Page 26
+                             _______________
06. DESCRIPTION OF MP SUBROUTINES
+________________________________
0We give first the method  of  calling  each  multiple-precision  routine
 directly,  second  (third,  ...)  alternative methods (if any) using the
 Augment interface described in Sections 4 and 5.
0Unless  otherwise  noted,  X, Y and Z denote integer arrays representing
 multiple-precision variables (often called (unpacked) multiple-precision
 numbers,  or  (unpacked)  multiple-precision  variables, or variables of
 type MULTIPLE), ERR and LV denote LOGICAL variables, I, J, K, L, IX etc.
 denote   single-precision   INTEGER   variables,   RX,  RY  etc.  denote
 single-precision REAL variables, DX, DY  etc.  denote  DOUBLE  PRECISION
 variables, and A denotes a packed Hollerith string.
0For definitions of B, T, M, LUN, MXR, RNDRL,  INBASE,  OUTBAS,  ...  see
 Section  3.1.   'Space'  means the size of the work area in COMMON, i.e.
 MXR+1-MNSPTR (see  Section  3.1).   MPGAM,  MPLNGM  and  MPZETA  require
 O(T**2) space, other MP routines require O(T) space.
0Time and space bounds such as O(T**2) are as T tends  to  infinity  with
 everything  else  fixed.   For  the  definition  of  the  function  m(T)
 appearing in some time bounds, see the description of MPMUL below.
0MPABS
+_____
0CALL MPABS (X, Y) or Y = ABS (X)
0Sets Y = abs(X) for multiple-precision variables X and Y.  X and  Y  may
 be  both  packed  or both unpacked.  This is an exception to the general
 rule that an unpacked result is returned even  if  the  argument(s)  are
 packed.   The result is exact.
0MPADD
+_____
0CALL MPADD (X, Y, Z) or Z = X + Y
0Adds   X   and   Y,   forming  result  in  Z,  where  X,  Y  and  Z  are
 multiple-precision numbers.
0Rounding  is  defined  by  the  parameter  RNDRL in COMMON /MPCOM/  (see
 Section 3.1) as follows -
0      RNDRL = 0 - Truncate towards zero if X*Y >= 0, away from  zero  if
                   X*Y  <  0, in both cases using one guard digit, so the
                   result is exact if severe cancellation occurs.
0      RNDRL = 1, 2  or  3  - See description of MPCQM.  Sufficient guard
                   digits are used to give the best possible result.
0MPADDI
+______
0CALL MPADDI (X, IY, Z) or Z = X + IY
0Adds  multiple-precision  X  to  integer IY giving multiple-precision Z.
 Rounding options are implemented as for MPADD.
1Sec. 6                       MP User's Guide                     Page 27
+                             _______________
0MPADDQ
+______
0CALL MPADDQ (X, I, J, Y) or Y = ADDQ (X, I, J)
0Adds the rational number I/J to the multiple-precision number X,  giving
 a  multiple-precision  result  in Y.  For the effect of rounding options
 see the descriptions of MPCQM and MPADD.
0MPADD3
+______
0Called  by  MPADD,  does  inner  loops of addition.  Not recommended for
 independent use.
0MPART1
+______
0CALL MPART1 (N, Y) or Y = ART1 (N)
0Computes multiple-precision Y = arctan(1/N), for integer N  >  1.   Uses
 the Taylor series arctan(x) = x - x**3/3 + x**5/5 - ... , where x = 1/N.
 Rounding options are implemented as for MPEXP.
0
 MPASIN
+______
0CALL MPASIN (X, Y) or Y = ASIN (X)
0Returns Y = arcsin(X), for multiple-precision numbers  X  and  Y,  where
 abs(X)  <=  1.  The result is in the range -pi/2 to +pi/2.  Uses MPATAN,
 so time = O(m(T)T).  Rounding options are not yet implemented, no  guard
 digits used.
0
 MPATAN
+______
0CALL MPATAN (X, Y) or Y = ATAN (X)
0Returns  Y  =  arctan(X)  for (packed or unpacked) multiple-precision X,
 unpacked  multiple-precision  Y,   using   an   O(m(T)T)   method.    An
 asymptotically  faster  method  would be to combine Newton's method with
 MPTAN.  For another asymptotically faster method, see  -  R.  P.  Brent,
 'Fast  multiple-precision evaluation of elementary functions', J. ACM 23
 (1976), 242-251, and the description of MPPIGL.  Result is in the  range
 -pi/2 to +pi/2.  Rounding options are implemented as for MPEXP.
0
 MPATN2
+______
0CALL MPATN2 (X, Y, Z) or Z = ATAN2 (X, Y)
0Returns Z = arctan (X/Y) if Y nonzero,
         Z = pi/2         if Y = 0,
0where  X  and  Y  are  (both packed or both unpacked) multiple-precision
 variables, Z  is  an  unpacked  multiple-precision  variable.   Rounding
 options are implemented as for MPEXP.
1Sec. 6                       MP User's Guide                     Page 28
+                             _______________
0MPBASA
+______
0I = MPBASA (X) or I = BASE (X)
0Returns  the multiple-precision base (first word of COMMON /MPCOM/ - see
 Section 3.1).  X is a dummy multiple-precision argument (Augment expects
 one), I is an integer.
0MPBASB
+______
0CALL MPBASB (I, X) or BASE (X) = I
0Sets  the  multiple-precision  base (first word of COMMON /MPCOM/) to I.
 I should be an integer such that I > 1 and (8*I*I-1) is representable as
 a   single-precision   integer   (see   Section  3.1).   X  is  a  dummy
 multiple-precision argument (Augment expects one).
0Setting  the  base  does  not  automatically  convert multiple-precision
 numbers to the new base.  This can be  done  by  converting  to  decimal
 using MPFOUT, changing the base, and converting back using MPFIN.
0MPBERN
+______
0CALL MPBERN (N, P, X) or X(1) = BERN (N, (PARAM ('MT2') + 1)/2)
0Computes  the  Bernoulli  numbers B(2) = 1/6, B(4) = -1/30, B(6) = 1/42,
 B(8)= -1/30, ... , B(2N), which are defined by the  generating  function
 y/(exp(y)-1).  N and P are single-precision integers, with 2*P > T+1.
0For use without Augment, X should be a one-dimensional integer array  of
 dimension  at  least P*N.  The Bernoulli numbers B(2), B(4), ... , B(2N)
 are returned in packed format in X, with B(2J) in locations
 X((J-1)*P+1),  ...  ,  X(J*P).   Thus,  to   get   B(2J)   in   unpacked
 multiple-precision  format in Y, one should CALL MPUNPK (X(IX), Y) after
 calling MPBERN, where IX = (J-1)*P+1.
0Alternatively  (simpler  but  nonstandard)  - X may be a two-dimensional
 integer array declared with dimension (P,  N1),  where   N1  >=  N   and
 2*P  >  T+1.  Then B(2), B(4), ... , B(2N) are returned in packed format
 in X, with B(2J) in X(1,J), ... , X(P,J).  Thus, to  get  B(2J)  in  the
 usual multiple-precision format in Y one should CALL MPUNPK (X(1, J), Y)
 after calling MPBERN.
0For use with Augment, declare
0       MULTIPAK X(N1)
0where N1 >= N, and use (PARAM('MT2') + 1)/2 as the second argument.
0The well-known recurrence is unstable (losing about 2J bits of  relative
 accuracy  in  the  computed  B(2J)),  so  we  use a different recurrence
 derived by multiplying the generating function y/(exp(y)-1) by sinh(y/2)
 and  equating  coefficients.   (This  method  was suggested by Christian
 Reinsch, and is faster than the  method  used  in  earlier  versions  of
 MPBERN.)
1Sec. 6                       MP User's Guide                     Page 29
+                             _______________
0The relation
0      B(2J) = -2((-1)**J)factorial(2J)zeta(2J)/((2pi)**(2J))
0is used if zeta(2J) equals 1 to working accuracy.  The relative error in
 B(2J) is O((J**2)*(B**(1-T))).  For further details, see - R. P.  Brent,
 'Unrestricted  algorithms  for  elementary   and   special   functions',
 Information Processing 80, North Holland, 1980, 613-619.
0Time = O(T*min(N,T)**2 + N*m(T)).
0Rounding options are not implemented, no guard digits used.
0The above remarks assume that N is positive.  If  N  is  negative,  they
 apply  with  N  replaced  by abs(N), except that the precision decreases
 linearly, from full precision for B(2) to low precision for B(2N).  This
 is  faster,  and  sufficiently  accurate for most applications where the
 Bernoulli  numbers  are  required  to  generate   coefficients   in   an
 Euler-Maclaurin sum formula.
0
0MPBESJ
+______
0CALL MPBESJ (X, NU, Y) or Y = BESJ (X, NU)
0Returns Y = J (NU, X), the first-kind Bessel function of order  NU,  for
 small  integer  NU,  multiple-precision  X  and  Y.   The method used is
 Hankel's asymptotic expansion if abs(X) is large, the  power  series  if
 abs(X)  is small, and the backward recurrence method otherwise.  Results
 for negative arguments are defined by
0      J (-NU, X) = J (NU, -X) = ((-1)**NU) * J (NU, X).
0Time = O(m(T)T) for fixed X and NU, increases  as  X  and  NU  increase,
 unless  X  is  large  enough  for  asymptotic  series to be used.  Space
 required is large (compared to most  multiple-precision  routines),  but
 O(T).
0Rounding options are not yet  implemented.   Uses  truncated  arithmetic
 with some guard digits.
0MPBES2
+______
0Uses the backward recurrence method to evaluate J (NU, X), where X is  a
 multiple-precision  variable, NU (the index) is an integer, and J is the
 Bessel function of the first kind.  Assumes that NU >= 0 and X > 0.  For
 normalization the identity
0      J(0,X) + 2*J(2,X) + 2*J(4,X) + ... = 1
0is  used.   Called  by  MPBESJ  and not recommended for independent use.
 Rounding options are not implemented, no guard digits used.
1Sec. 6                       MP User's Guide                     Page 30
+                             _______________
0MPCAM
+_____
0CALL MPCAM (A, X) or X = CTM (A) or X = CAM (A)
0Converts the  (packed) Hollerith string A to a multiple-precision number
 X.   A  can  be  a  string  of  digits  acceptable  to  routine MPIN and
 terminated by '$', e.g. '-5.367$' or '1E-99$', or one of  the  following
 special strings
0      'EPS'   (multiple-precision machine-precision, see MPEPS),
       'EUL'   (Euler's constant 0.5772..., see MPEUL),
       'MAXR'  (largest valid multiple-precision number, see MPMAXR),
       'MINR'  (smallest positive multiple-precision number, see MPMINR),
       'PI'    (pi = 3.14..., see MPPI).
0Actually,  only  the  first  two characters of these special strings are
 significant. For example, 'MA', 'MAX' and 'MAXR' are equivalent.
0The  error  message  '***  MXR  TOO  SMALL ... ***' is usually caused by
 omission of the sentinel '$'.
0For the effect or rounding options, see the descriptions of MPIN, MPEPS,
 MPMAXR, MPMINR and MPPI.
0Warning to Augment users - use CAM(A) and not CTM(A) if A is declared as
 an integer array.  (Otherwise Augment will  generate  a  call  to  MPCIM
 instead of MPCAM.)
0MPCDM
+_____
0CALL MPCDM (DX, Z) or Z = DX or Z = CTM (DX)
0Converts double-precision DX to multiple-precision Z.  Some numbers will
 not convert exactly on machines  with  base  other  than  two,  four  or
 sixteen,  or  if  B  is not a power of two, or if T is too small.  Thus,
 MPCDM should be used only to obtain starting  approximations  etc.,  and
 not  where  full  multiple-precision  accuracy  is  required.   Rounding
 options  are  not   implemented.    For   accurate   initialization   of
 multiple-precision variables, use MPCAM, MPCIM, MPCQM or MPQPWR.
0
 MPCDM, MPCMD, MPCMDE and MPCMPD are not called by any other routines  in
 the  MP  package,  so  they  may  be  omitted if double-precision is not
 available.
0
 MPCEIL
+______
0CALL MPCEIL (X, Y) or Y = CEIL (X)
0Sets Y = ceiling (X), i.e. the smallest integer not less than X, where X
 and Y are unpacked multiple-precision variables.  Rounding is defined by
 RNDRL in COMMON /MPCOM/, as for MPEXP (this is only  relevant  if  X  is
 large  and  positive, for otherwise ceiling (X) is exactly representable
 as a multiple-precision number).
1Sec. 6                       MP User's Guide                     Page 31
+                             _______________
0MPCHEB
+______
0CALL MPCHEB (C, NC, N, IND) or CALL MPCHEB (C, PARAM ('MT2'), N, IND)
0Converts the power series coefficients C(1), ...  ,  C(N)  to  Chebyshev
 series  coefficients  (which overwrite C(1), ... , C(N)).  It is assumed
 that on summation of the Chebyshev series the constant term is halved.
0IND may have the value 0, -1 or +1, with the following meanings -
0  0 - C(j)  is the coefficient of x**(j-1) on input, and of T(j-1, x) on
       output, where T(0, x), T(1, x), ... are the Chebyshev polynomials,
       e.g.
             T(3, x) = 4*x**3 - 3*x.
0 -1 - C(j)  is  the coefficient of  x**(2j-1)  on input,  and of T(2j-1,
       x) on output (useful for approximation of odd functions).
0 +1 - C(j)  is  the coefficient of  x**(2j-2)  on input,  and of T(2j-2,
       x) on output (useful for approximation of even functions).
0Conceptually, C is an array of multiple-precision variables.  If  MPCHEB
 is  called  directly,  C  is a two-dimensional INTEGER array (with first
 dimension NC at least T+2, second dimension at least N).  If Augment  is
 used,  C  should  be of type MULTIPLE, with dimension at least N, and NC
 should be PARAM ('MT2').
0Rounding  options are not implemented.  Uses truncated arithmetic and no
 guard digits, time = O(N*N*T), space = O(T).
0MPCHEV
+______
0CALL MPCHEV (C, NC, N, IND, X, Y)
0                         or CALL MPCHEV (C, PARAM ('MT2'), N, IND, X, Y)
0Returns  Y = Chebyshev series evaluated at X, where X and Y are unpacked
 multiple-precision variables, C is an  array  containing  the  Chebyshev
 coefficients.   For a description of the arguments C, NC, N and IND, see
 the notes on MPCHEB above.
0The  algorithm  used  is  described  in  Chapter 8 of - Modern Computing
 Methods (Notes on Applied Science No. 16, HMSO, London, 1961).
 TIme = O(N*m(T)), space = O(T).
0Rounding  options are not implemented.  Uses truncated arithmetic and no
 guard digits.
0MPCHGB
+______
0J = MPCHGB (B1, B2, N)
0Returns J such that B1**abs(J) >= B2**abs(N), i.e.
0      abs(J) >= abs(N)*log(B2)/log(B1),
1Sec. 6                       MP User's Guide                     Page 32
+                             _______________
0and sign(J) = sign(N), where B1, B2 and J are integers, B1 > 1, B2 >  0,
 B1*B2 <= MXINT.  The value of J returned is close to minimal.  Uses only
 integer arithmetic.
0MPCHK
+_____
0CALL MPCHK
0Checks legality of parameters in COMMON /MPCOM/ (see Section  3.1),  and
 updates  MXSPTR.  If an illegal parameter in COMMON /MPCOM/ is detected,
 an error message is written on unit LUN, and execution is terminated.
0MPCIM
+_____
0CALL MPCIM (I, Z) or Z = I or Z = CTM (I)
0Converts single-precision  integer  I  to  multiple-precision  Z,  where
 abs(I) <= B**T.
0MPCIS
+_____
0CALL MPCIS (X, C, S, BOTH)
0Returns  C = cos (X) and S = sin (X)  if BOTH = .TRUE.,
 or just  C = cos (X)                  if BOTH = .FALSE.
0X, C and S are unpacked multiple-precision variables, BOTH is a  logical
 variable.   Faster  than  calling  MPCOS  and  MPSIN with the same first
 argument.  The  algorithm  is  similar  to  that  used  in  MPEXP  (with
 imaginary  argument),  time  =  O(sqrt(T)m(T)).   Rounding  options  are
 implemented as follows -
0      RNDRL = 0 or 1 - Absolute error < 0.6*B**(-T)
                     (but relative error may be large if abs(X) > 1).
       RNDRL = 2  -  Lower bound on the true result.
       RNDRL = 3  -  Upper bound on the true result.
0MPCMD
+_____
0CALL MPCMD (X, DZ) or DZ = X or DZ = CTD (X)
0Converts multiple-precision X to double-precision DZ.  Assumes X  is  in
 allowable range for double-precision numbers (otherwise double-precision
 floating-point overflow or underflow may occur).  There may be some loss
 of  accuracy  if  B is not a power of the base used for double-precision
 floating-point arithmetic, or if T is too small.  Rounding  options  are
 not implemented.  See also the description of MPCDM.
0MPCMDE
+______
0CALL MPCMDE (X, N, DX)
0Returns integer N and double-precision DX such that (multiple-precision)
0      X = DX*OUTBAS**N (approximately),
1Sec. 6                       MP User's Guide                     Page 33
+                             _______________
0where 1 <= abs(DX) < OUTBAS (unless X = 0).  The default value of OUTBAS
 is ten (see Section 3.1).  It is assumed that X is not so large or small
 that  N  overflows.   X may be packed or unpacked.  Rounding options are
 not implemented.  See also the description of MPCDM.
0
 MPCMEF
+______
0CALL MPCMEF (X, N, Y)
0Given multiple-precision X, returns integer N and  multiple-precision  Y
 such that
0      X = (OUTBAS**N)*Y and 1 <= abs(Y) < OUTBAS
0(unless X = 0, when N = Y = 0).  The default value of OUTBAS is ten (see
 Section 3.1).  It is assumed that X is not so  large  or  small  that  N
 overflows.
0Rounding options are not fully implemented, but  the  directed  rounding
 options  (RNDRL  = 2 or 3) give correct (though not best possible) lower
 and upper bounds on the true result.
0
 MPCMF
+_____
0CALL MPCMF (X, Y) or Y = CMF (X) or Y = FRAC (X)
0For multiple-precision X and Y, returns fractional part of X in Y,  i.e.
0      Y = X - integer part of X (truncated towards 0).
0Rounding options are irrelevant as the result is exact.
0
 MPCMI
+_____
0CALL MPCMI (X, I) or I = X or I = CTI (X)
0Converts multiple-precision X to integer I, assuming that X is  not  too
 large  (else  use MPCMIM).  X is truncated towards zero (irrespective of
 RNDRL).  I is returned as zero if abs(int(X)) >  MXINT.   The  user  may
 check for this possibility by testing if
0      ((X(1).NE.0) .AND. (X(2).GT.0) .AND. (I.EQ.0))
0is true on return from MPCMI.
0MPCMIM
+______
0CALL MPCMIM (X, Y) or Y = CMIM (X) or Y = INT (X)
0Returns   Y   =   integer   part   of   X  (truncated  towards  0),  for
 multiple-precision X and Y.  Use if abs(Y) > MXINT (otherwise  MPCMI  is
 preferable.)   X may be packed or unpacked if MPCMIM is called directly.
1Sec. 6                       MP User's Guide                     Page 34
+                             _______________
0MPCMP
+_____
0J = MPCMP (X, Y)
0Compares the unpacked multiple-precision numbers X and Y, returning
0      +1 if X > Y,
        0 if X = Y,
       -1 if X < Y.
0X  and  Y  must  be unpacked.  If they are packed, MPCOMP should be used
 instead.  The result is exact.
0
 MPCMPA
+______
0J = MPCMPA (X, Y) or J = CMPA (X, Y)
0Compares abs(X) with abs(Y) for multiple-precision X and Y, returning
0      +1 if abs(X) > abs(Y),
        0 if abs(X) = abs(Y),
       -1 if abs(X) < abs(Y).
0X and/or Y may be packed or unpacked if MPCMPA is called  directly,  but
 must be unpacked if called via Augment.  The result is exact.
0
 MPCMPD
+______
0J = MPCMPD (X, DI) or J = COMP (X, DI)
0Compares multiple-precision X with double-precision DI, returning
0      +1 if X > DI,
        0 if X = DI,
       -1 if X < DI.
0X may be packed or unpacked if MPCMPD is called directly.  Uses MPCDM to
 convert DI to multiple-precision, so the  result  may  be  incorrect  if
 MPCDM is inaccurate.  See the description of MPCDM.
0
 MPCMPI
+______
0J = MPCMPI (X, I) or J = COMP (X, I)
0Compares multiple-precision number X with integer I, returning
0      +1 if X > I,
        0 if X = I,
       -1 if X < I.
0X  may  be  packed or unpacked if MPCMPI is called directly, but must be
 unpacked if called via Augment.  The result is exact.
1Sec. 6                       MP User's Guide                     Page 35
+                             _______________
0MPCMPQ
+______
0K = MPCMPQ (X, I, J) or K = COMP (X, I, J)
0Compares multiple-precision X with the rational number I/J, where I  and
 J are single-precision integers, J nonzero.  Returns
0      +1 if X > I/J,
        0 if X = I/J,
       -1 if X < I/J.
0Rounding  options  are  irrelevant  as the result is exact (the rational
 number I/J is never computed).
0
0MPCMPR
+______
0J = MPCMPR (X, RI) or J = COMP (X, RI)
0Compares multiple-precision number X with real number RI, returning
0      +1 if X > RI,
        0 if X = RI,
       -1 if X < RI.
0X may be packed or unpacked if MPCMPR is called directly.  Uses MPCRM to
 convert RI to multiple-precision, so the  result  may  be  incorrect  if
 MPCRM is inaccurate.  See the description of MPCRM.
0
0MPCMR
+_____
0CALL MPCMR (X, RZ) or RZ = X or RZ = CTR (X)
0Converts multiple-precision X to single-precision real RZ.  Assumes that
 X  is  in  the  allowable  range  for  single-precision  real  variables
 (otherwise floating-point overflow or underflow may occur).   There  may
 be  some  loss  of  accuracy  if  B  is not a power of the base used for
 single-precision real arithmetic.  Rounding options are not implemented.
0
 MPCMRE
+______
0CALL MPCMRE (X, N, RX)
0Returns   integer   N   and   single-precision   real   RX   such   that
 multiple-precision
0      X = RX*OUTBAS**N (approximately),
0where 1 <= abs(RX) < OUTBAS (unless X = 0).  It is assumed that X not so
 large or small that N overflows.  The default value  of  OUTBAS  is  ten
 (see Section 3.1).  Rounding options are not implemented.
1Sec. 6                       MP User's Guide                     Page 36
+                             _______________
0MPCMUL
+______
0CALL MPCMUL (XR, XI, YR, YI, ZR, ZI)
0Computes the complex product
0      ZR + i*ZI = (XR + i*XI) * (YR + i*YI),
0where i**2 = -1.  XR, XI, YR, YI, ZR, ZI are unpacked multiple-precision
 variables.   Rounding  options  are  not  implemented.   Uses  truncated
 arithmetic with no guard digits.  Time = O(m(T)).
0MPCOMP
+______
0J = MPCOMP (X, Y) or J = COMP (X, Y)
0Compares the multiple-precision numbers X and Y, returning
0      +1 if X > Y,
        0 if X = Y,
       -1 if X < Y.
0X  and/or  Y may be packed or unpacked if MPCOMP is called directly, but
 must be unpacked if called via Augment.  The result is exact.
0MPCOS
+_____
0CALL MPCOS (X, Y) or Y = COS (X)
0Returns Y = cos(X) for multiple-precision X and Y, using MPCIS.   X  may
 be  packed or unpacked, Y is unpacked.  Time = O(sqrt(T)m(T)).  Rounding
 options are implemented as for MPCIS, so the absolute error is small but
 the relative error may be large.
0MPCOSH
+______
0CALL MPCOSH (X, Y) or Y = COSH (X)
0Returns multiple-precision Y = COSH(X) for multiple-precision X (not too
 large).  X may be packed or unpacked, Y is unpacked.   Rounding  options
 are not yet implemented, uses no guard digits.
 Time = O(sqrt(T)m(T)).
0MPCQM
+_____
0CALL MPCQM (I, J, Q) or Q = CQM (I, J) or Q = CTM (I, J)
0Converts  the  rational  number  I/J  to multiple-precision Q.  Rounding
 options are implemented.  Rounding is controlled by the parameter  RNDRL
 in COMMON /MPCOM/ (see Section 3.1), as follows -
0  0 - Round towards zero (i.e. chop or truncate). The resulting error is
       less than 1 ulp (unit in the last place).   This  is  the  default
       option if MPINIT or MPSET2 is used.  It is the fastest option, and
       is satisfactory for most purposes.
1Sec. 6                       MP User's Guide                     Page 37
+                             _______________
0
   1 - Round to nearest, preferring even last  digit  if  a  tie  occurs.
       Error  is  less  than  or  equal  to  0.5  ulp.  This option gives
       worst-case error bounds half as  large  as  for  RNDRL  =  0,  and
       average-case  error  considerably  better  than  for  RNDRL  =  0.
       However, it is expensive in time and space.  If RNDRL = 0 does not
       give  sufficient accuracy, it is usually better to increase T than
       to use RNDRL = 1.
0  2 - Round  down  (towards -infinity).  Error is less than 1 ulp.  This
       option is useful for interval arithmetic.
0  3 - Round  up  (towards  +infinity).   Error is less than 1 ulp.  This
       option is useful for interval arithmetic.
0MPCRM
+_____
0CALL MPCRM (RX, Z) or Z = RX or Z = CTM (RX)
0Converts single-precision real RX to multiple-precision Z.  Some numbers
 will not convert exactly on machines with base other than two,  four  or
 sixteen, or if B is not a power of 2, or if T is too small.  Thus, MPCRM
 should only be used to generate starting approximations  etc.,  and  not
 where  full  multiple-precision  accuracy  is  required.   For  accurate
 initialization of multiple-precision variables, use MPCAM, MPCIM,  MPCQM
 or MPQPWR.  Rounding options are not implemented.
0MPDAW
+_____
0CALL MPDAW (X, Y) or Y = DAW (X)
0Returns
0      Y = Dawson's integral of multiple-precision argument X
0        = exp(-X**2)*(integral from 0 to X of exp(u**2)du),
0where X and Y are multiple-precision variables.   X  may  be  packed  or
 unpacked  if MPDAW is called directly.  Rounding options are implemented
 as for MPEXP.
0MPDGA
+_____
0I = MPDGA (X, N) or I = DIGIT (X, N)
0Returns the N-th digit of the multiple-precision number X,  0 < N <=  T.
 Returns zero if X is zero.
0MPDGB
+_____
0CALL MPDGB (I, X, N) or DIGIT (X, N) = I
0Sets  the N-th digit of the multiple-precision number X to I.  N must be
 in the range 0 < N <= T, I must be in the range 0 <= I < B (and I > 0 if
 N = 1).  The sign and exponent of X are not changed.
1Sec. 6                       MP User's Guide                     Page 38
+                             _______________
0MPDIGA
+______
0I = MPDIGA (X) or I = NUMDIG (X)
0Returns  T,  the number of multiple-precision digits (the second word of
 COMMON /MPCOM/).  X is  a  dummy  multiple-precision  argument  (Augment
 expects one).
0MPDIGB
+______
0CALL MPDIGB (I, X) or NUMDIG (X) = I
0Sets  T,  the number of multiple-precision digits (second word of COMMON
 /MPCOM/), to I, which should be an integer in the range 1 < I  <  MT2-1.
 X is a dummy multiple-precision argument (Augment expects one).
0MPDIGS
+______
0J = MPDIGS (N)
0Returns  the  number  of multiple-precision (base B) digits required for
 the equivalent of at least N floating 'decimal' places.   The inequality
0      B**(MPDIGS - 1) >= OUTBAS**(N - 1)
0will be satisfied, usually with the minimal such value  of  MPDIGS.   If
 OUTBAS  has its default value (ten), the working precision may be set to
 the equivalent of at least N floating decimal places by the statement
0      PARAM ('NUMDIG') = MPDIGS (N)       (using Augment)
 or
       CALL MPPARB (MPDIGS (N), 'NUMDIG')  (not using Augment).
0MPDIGV
+______
0J = MPDIGV (C)
0Returns the integer value
0      0 ,  ... ,  9 ,  10, ... ,  15 if C is the corresponding character
      '0',  ... , '9', 'A', ... , 'F'
0and  the  value  returned  is  less  than INBASE (default value ten, see
 Section 3.1).  Otherwise returns -1.   The  character  value  should  be
 represented in an integer variable as if read under A1 format.
0MPDIGW
+______
0C = MPDIGW (J)
0Returns the character '0', ... , 'F'  if J is the corresponding integer,
                        0 , ... , 15,
0returns '*' otherwise.  The character is stored in an integer as if read
 under A1 format.
1Sec. 6                       MP User's Guide                     Page 39
+                             _______________
0MPDIM
+_____
0CALL MPDIM (X, Y, Z) or Z = DIM (X, Y)
0Sets Z = X - min (X, Y) = max (0, X-Y) for multiple-precision  variables
 X, Y and Z.  X and Y may be both packed or both unpacked, Z is unpacked.
 Rounding options are implemented as for MPSUB.
0
 MPDIV
+_____
0CALL MPDIV (X, Y, Z) or Z = X/Y
0Sets Z = X/Y, for multiple-precision X, Y and Z.  Y must  not  be  zero.
 Uses  MPDIVL  for  small  T  or if RNDRL > 0, MPREC and MPMUL otherwise.
 Thus time = O(m(T) + RNDRL*T**2).  Rounding options are  implemented  as
 for MPCQM.
0
 MPDIVI
+______
0CALL MPDIVI (X, IY, Z) or Z = X/IY
0Divides  unpacked  multiple-precision  X by the nonzero single-precision
 integer IY, giving multiple-precision Z in  time  O(T).   This  is  much
 faster  than  division by a multiple-precision number.  Rounding options
 are implemented as for MPCQM.
0MPDIVL
+______
0CALL MPDIVL (X, Y, Z)
0Sets  Z = X/Y  for  unpacked  multiple-precision  variables  X,  Y and Z
 (Y  nonzero).     Uses  the 'schoolboy' or 'long division' algorithm, so
 time = O(T**2).  An alternative method is to use MPREC and MPMUL.   This
 would be faster for large T, but MPDIVL is faster for small T.  Rounding
 options are implemented as for MPCQM.
0MPDIV2
+______
0Routine called by MPDIVI, not recommended for independent use.
0MPDIV3
+______
0Routine called by MPDIVL, not recommended for independent use.
0MPDUMP
+______
0CALL MPDUMP (X)
0Dumps the (packed or unpacked) multiple-precision number X on unit  LUN.
 Writes  the  sign, exponent and digits in suitable I format (or just the
 sign if X is zero).  Embedded blanks should  be  interpreted  as  zeros.
 MPDUMP  is  sometimes  useful  for  debugging, but is not recommended in
 general for output of MP numbers.
1Sec. 6                       MP User's Guide                     Page 40
+                             _______________
0MPEI
+____
0CALL MPEI (X, Y) or Y = EI (X)
0Returns
0      Y = ei (X) = -E1 (X)
         = (principal value integral from -infinity to X of exp(u)/u du),
0for multiple-precision numbers X and Y, using the power series for small
 abs(X), the asymptotic series for  large  abs  (X),  and  the  continued
 fraction for medium negative X.  The relative error in Y is small unless
 X is very close to the zero 0.37250741078136663446... of ei(X), and then
 the  absolute  error  in  Y  is O(B**(1-T)).  In any case the error in Y
 could be induced by an O(B**(1-T)) relative perturbation in X.  Rounding
 options are not yet implemented.  Time = O(m(T)T).
0
0MPEPS
+_____
0CALL MPEPS (X) or X = CTM ('EPS')
0Sets multiple-precision X to the (multiple-precision) machine precision,
 that  is  an upper bound on the smallest representable positive number X
 such that the relative error in the basic multiple-precision  operations
 of  addition, subtraction, multiplication and division does not exceed X
 unless the result underflows.  In fact
0      X = 1.01*B**(1-T)  (rounded up)  if RNDRL = 0,
         =  0.5*B**(1-T)  (rounded up)  if RNDRL = 1,
         =      B**(1-T)                if RNDRL = 2 or 3.
0
0MPEQ
+____
0LV = MPEQ (X, Y) or LV = (X .EQ. Y) or IF (X .EQ. Y) ...
0Returns logical value LV of (X .EQ. Y) for multiple-precision X  and  Y.
 MPEQ  must  be  declared  LOGICAL  unless the Augment interface is used.
 X and/or Y may be packed or unpacked.
0
0MPERF
+_____
0CALL MPERF (X, Y) or Y = ERF (X)
0Returns Y = erf(X)
           = sqrt(4/pi)*(integral from 0 to X of exp(-u**2) du)
0for multiple-precision X and Y.  X may be packed or unpacked if MPERF is
 called directly, must be  unpacked  if  MPERF  is  called  via  Augment.
 Rounding options are implemented as for MPEXP.
1Sec. 6                       MP User's Guide                     Page 41
+                             _______________
0MPERFC
+______
0CALL MPERFC (X, Y) or Y = ERFC (X)
0Returns Y = erfc(X) = 1 - erf(X) for multiple-precision numbers X and Y.
 X may be packed or unpacked  if  MPERFC  is  called  directly,  must  be
 unpacked if MPERFC is called via Augment.  Rounding as for MPEXP.
0MPERF2
+______
0Computes exp(X**2)*(integral from 0 to X of exp(-u*u) du)
0for  multiple-precision X, using the power series for small X, and MPEXP
 for large X.  Called by MPERF,  not  recommended  for  independent  use.
 Rounding options are not implemented, no guard digits are used.
0MPERF3
+______
0Trys to return
0      exp(X**2)*(integral from X to infinity of exp(-u**2)du),
 or
       exp(-X**2)*(integral from 0 to X of exp(u**2) du),
0using   the  asymptotic  expansion,  where  X  is  a  multiple-precision
 variable.  The condition for success is approximately that
0      X > sqrt(T*ln(B)).
0Called by MPERF, MPERFC and MPDAW, not recommended for independent  use.
 Rounding options are not implemented, uses no guard digits.
0MPERR
+_____
0CALL MPERR
0This  routine is called when a fatal error condition is encountered, and
 after a message has been written on logical unit LUN.   As  supplied  in
 the  multiple-precision  package,  MPERR writes a message which includes
 the version number in the format YYMMDD, then stops.  On many systems it
 can  easily  be  modified  to  give  a trace-back (e.g. with Univac 1100
 Fortran V this can be obtained by replacing STOP by RETURN 0).
0Since  MPERR  is  only  called  when  a  fatal  error condition has been
 detected, it does not return to the calling routine.
0MPERRM
+______
0CALL MPERRM (A)
0The argument A is a packed Hollerith string terminated by  '$'.   MPERRM
 writes  the  string  on  unit  LUN, then calls MPERR.  Only the first 71
 characters of A are significant.  They will be preceeded and followed by
 '***'.   MPERRM  is  called  by  MP  routines  when  a  fatal  error  is
 encountered, and does not return to the calling routine.
1Sec. 6                       MP User's Guide                     Page 42
+                             _______________
0MPEUL
+_____
0CALL MPEUL (G) or G = CTM ('EUL')
0Returns multiple-precision G = Euler's constant (gamma =  0.57721566...)
 to  almost  full multiple-precision accuracy.  The method was discovered
 by Edwin McMillan and Richard Brent, and is faster than  the  method  of
 Sweeney  (used  in earlier versions of MPEUL).  See - R. P. Brent and E.
 M. McMillan, 'Some new  algorithms  for  high-precision  computation  of
 Euler's  constant',  Math.  Comp.  34(1980),  305-312.   Time = O(T**2).
 Rounding options are implemented as for MPEXP.
0MPEXP
+_____
0CALL MPEXP (X, Y) or Y = EXP (X)
0Returns Y = exp(X) for multiple-precision X and Y.  X may be  packed  or
 unpacked if MPEXP is called directly.  The method used is described in -
 R. P. Brent, 'The complexity of multiple-precision arithmetic',  in  The
 Complexity  of  Computational  Problem  Solving,  Queensland  University
 Press, Brisbane, 1976, 126-165.  See also - R. P.  Brent,  'Unrestricted
 algorithms for elementary and special functions', Information Processing
 80, North-Holland, 1980, 613-619.  Time = O(sqrt(T)m(T)).
0Rounding  is  controlled  by  the parameter RNDRL in COMMON /MPCOM/ (see
 Section 3.1), as follows -
0      0 or 1 - error less than 0.6 ulp (units in the last place),
0      2      - computed result is lower bound on the true result,
0      3      - computed result is upper bound on the true result.
0
 MPEXPA
+______
0I = MPEXPA (X) or I = EXPON (X)
0Returns  the  exponent  of  the  multiple-precision  number  X (or large
 negative exponent -M if X is zero).  X may  be  packed  or  unpacked  if
 MPEXPA is called directly.
0MPEXPB
+______
0CALL MPEXPB (I, X) or EXPON (X) = I
0Sets  the  exponent of the multiple-precision number X to I, unless X is
 zero  (when  its  exponent  is  unchanged).    X   must   be   a   valid
 multiple-precision  number (either zero or normalized).  X may be packed
 or unpacked if MPEXPB is called directly.  Note that  multiple-precision
 underflow/overflow will occur if I is too small/large.
0MPEXP1
+______
0CALL MPEXP1 (X, Y) or Y = EXP1 (X)
1Sec. 6                       MP User's Guide                     Page 43
+                             _______________
0Returns
0      Y = exp(X) - 1
0where  X  and  Y are multiple-precision numbers and -1 < X < 1.  Uses an
 O(sqrt(T)m(T) algorithm described in - R. P. Brent, 'The  complexity  of
 multiple-precision  arithmetic',  in Complexity of Computational Problem
 Solving, Univ. of Queensland Press, Brisbane, 1976,  126-165.   Rounding
 options are implemented as for MPEXP.
0MPFIN
+_____
0CALL MPFIN (LUNIT, X) or X = FIO (LUNIT)
0Reads  (unpacked)  multiple-precision number X from Fortran logical unit
 LUNIT  (LUNIT >= 0).  It is assumed that a record of length INRECL <= 80
 may  be  read  in '(80A1)' format from unit LUNIT, and then converted to
 multiple-precision  using  MPIN.   Thus,  MPFIN  may  be  used  to  read
 free-format  fixed  or floating-point numbers, one per record, from unit
 LUNIT.  The input base INBASE has default value ten (i.e. decimal  input
 - see Section 3.1).  For further details see under MPIN.
0MPFLOR
+______
0CALL MPFLOR (X, Y) or Y = FLOOR (X)
0Sets  Y  =  floor (X), i.e. the largest integer not exceeding X, where X
 and Y are unpacked multiple-precision variables.  Rounding is defined by
 RNDRL  in  COMMON  /MPCOM/,  as for MPEXP (this is only relevant if X is
 large and negative, for otherwise floor (X) is exactly representable  as
 a multiple-precision number).
0
 MPFOUT
+______
0CALL MPFOUT (X, N) or FIO (N) = X
0Writes  the  multiple-precision  number  X on unit LUN in floating-point
 format, with N significant 'decimal' digits (N > 1).  The exponent field
 width  is  EXWID  (default  is 6, including exponent character EXPCH and
 sign - see Section 3.1).  The default exponent character  EXPCH  is  'E'
 (or '$' if OUTBAS = 15 or 16), and the default output base OUTBAS is ten
 (i.e. decimal output).  Rounding options are implemented as for  MPOUTE.
 X  may  be  packed or unpacked if MPFOUT is called directly, but must be
 unpacked if called via Augment.
0MPGAM
+_____
0CALL MPGAM (X, Y) or Y = GAM (X)
0Computes multiple-precision Y = Gamma(X) for multiple-precision argument
 X,  using  MPGAMQ  if  240*X is a small integer, otherwise using MPLNGM.
 Space  required  is  about  the  same  as  for  MPLNGM,  i.e.   O(T**2),
 time  =  O(T**3).  Rounding options are not yet implemented and no guard
 digits are used.
1Sec. 6                       MP User's Guide                     Page 44
+                             _______________
0MPGAMQ
+______
0CALL MPGAMQ (I, J, X) or Y = GAMQ (I, J) or Y = GAM (I, J)
0Returns
0      X = Gamma (I/J),
0where X is multiple-precision, I and J are small integers.   The  method
 used  is reduction of the argument to (0, 1) and then a direct expansion
 of the defining integral truncated at a sufficiently high  limit,  using
 2T  digits to compensate for cancellation.  Time = O(T**2) if I/J is not
 too large.  If I/J > 100 (approximately) it  is  faster  to  use  MPGAM.
 (MPGAMQ  is  very  slow  if  I/J  is is very large, because the relation
 Gamma(X+1) = X*Gamma(X) is used repeatedly.)  Rounding options  are  not
 yet implemented.
0MPGCD
+_____
0CALL MPGCD (K, L)
0Returns  K  =  K/GCD  and  L  =  L/GCD, where GCD is the greatest common
 divisor of the initial K and L (single-precision integers).   Called  by
 various MP routines.
0
 MPGCDA
+______
0CALL MPGCDA (X, Y, Z) or Z = GCD (X, Y)
0Returns Z = greatest common divisor of X and Y.  (By definition,
 GCD  (X,  0)  =  GCD  (0, X) = abs(X), GCD (X, Y) >= 0.)  X, Y and Z are
 integers represented as multiple-precision  numbers,  and  must  satisfy
 abs(X)  <  B**T,  abs(Y)  <  B**T.   The  method  is  a straight-forward
 implementation of the Euclidean algorithm, and  Lehmer's  trick  is  not
 used,  although  it might be faster (see Amer. Math.  Monthly 45 (1938),
 227-233).  Time = O(T**2).  X and/or Y may  be  packed  or  unpacked  if
 MPGCDA  is  called  directly.   Rounding  options  are irrelevant as the
 result is exact.
0MPGCDB
+______
0CALL MPGCDB (X, Y)
0Returns (X, Y) as (X/Z, Y/Z) where Z is the GCD  of  initial  X  and  Y,
 which  are  integers represented as multiple-precision numbers, and must
 satisfy abs(X) < B**T, abs(Y) < B**T.  Time = O(T**2).  Rounding options
 are irrelevant as the result is exact.
0
 MPGD
+____
0J = MPGD (N)
0Returns ceiling (ln (max (1, abs(N))) / ln(B)),
1Sec. 6                       MP User's Guide                     Page 45
+                             _______________
0i.e. the minimum J >= 0 such that
0      B**J >= abs(N).
0This function is  useful  for  computing  the  number  of  guard  digits
 required for various (multiple-precision) base B calculations.
0MPGD3
+_____
0CALL MPGD3 (N, TG)
0Sets  T = T + 1 + MPGD (100*N)  if it is safe to call MPGD, and the same
 or slightly more otherwise.  Also sets TG to the new value of T.  N  and
 TG  are  integer  arguments.   Called  by various MP routines when it is
 necessary to increase the working precision.
0MPGE
+____
0LV = MPGE (X, Y) or LV = (X .GE. Y) or IF (X .GE. Y) ...
0Returns the logical value of (X > Y) for  multiple-precision  X  and  Y.
 MPGE  must  be  declared  LOGICAL  unless the Augment interface is used.
 X and/or Y may be packed or unpacked.
0MPGET
+_____
0J = MPGET (X, N)
0Returns X(N), where X is an integer  array  of  dimension  at  least  N.
 Necessary to avoid some compiler diagnostics.
0MPGT
+____
0LV = MPGT (X, Y) or LV = (X .GT. Y) or IF (X .GT. Y) ...
0Returns  the  logical  value  of (X > Y) for multiple-precision X and Y.
 MPGT must be declared LOGICAL unless the Augment interface is  used.   X
 and/or Y may be packed or unpacked.
0MPHANK
+______
0Tries to compute the Bessel function J (NU, X) using Hankel's asymptotic
 expansion,  for nonnegative integer NU and multiple-precision X.  Time =
 O(T**3).   Called  by  MPBESJ,  not  recommended  for  independent  use.
 Rounding options are not implemented, uses no guard digits.
0MPIMUL
+______
0CALL MPIMUL (IY, X, Z) or Z = IY*X
0Equivalent to
0      CALL MPMULI (X, IY, Z) or Z = X*IY.
0See the description of MPMULI for further details.
1Sec. 6                       MP User's Guide                     Page 46
+                             _______________
0MPIN
+____
0CALL MPIN (C, X, N, ERROR)
0Converts  the characters (assumed to have been read under NA1 format) in
 C(1), ... , C(N) to a multiple-precision number in X.  If C represents a
 valid  number,  ERROR  is returned as .FALSE.  If C does not represent a
 valid number, ERROR is returned as .TRUE. and X as  zero.   Leading  and
 trailing  blanks are allowed, embedded blanks (except between the number
 and its sign) are forbidden.  If  there  is  no  decimal  point  one  is
 assumed  to  lie  just  to  the right of the last decimal digit.  X is a
 multiple-precision number, C an integer array, N an integer,  and  ERROR
 logical.
0The input base is determined by the parameter INBASE in  COMMON  /MPCOM/
 (see  Section 3.1).  INBASE has default value ten, but may be set to any
 value in two, ... , sixteen.
0    Examples of valid numbers    Examples of invalid numbers
+    _________________________    ___________________________
0         -   123456789                    12 345
              3.14159                      123.456E -67
             -44.                          1.2.3
             .0001234                      E123
              123.456D789                  64.4E+
            -.1234566-789                ++12.3
             +999+88                       E3.
0Rounding is determined by the  parameter  RNDRL  in  COMMON  /MPCOM/  as
 follows -
0      0 or 1 - Round  to  (approximately) the nearest representable base
                B, T-digit number.  The error is less than 0.6  units  in
                the last (base B) place.
0      2   -    Round down (towards  -infinity)  to  a  base  B,  T-digit
                number,
0      3   -    Round up (towards +infinity) to a base B, T-digit number.
0MPINE
+_____
0CALL MPINE (C, X, N, J, ERROR)
0MPINE  is  the  same as MPIN except that the result (X) is multiplied by
 INBASE**J, where J is a single-precision integer.  For  details  of  the
 other   arguments,   see  MPIN.   Useful  for  floating-point  input  of
 multiple-precision numbers.  The user  can  read  the  exponent  into  J
 (using  any suitable format) and the fraction into C (using NA1 format),
 then call MPINE to convert to multiple-precision.
0Note - in early versions of the MP package, MPIN was unable to deal with
 numbers with exponents (e.g. 5.3E-26).  Now that MPIN can deal with such
 numbers  there  is  no  real  need for MPINE, so it is included only for
 compatibility with early versions of the package.
1Sec. 6                       MP User's Guide                     Page 47
+                             _______________
0MPINF
+_____
0CALL MPINF (X, N, UNIT, IFORM, ERR)  or  ERR = MPINF (X, N, UNIT, IFORM)
0                                   or IF (MPINF (X, N, UNIT, IFORM)) ...
0Reads  N  words  from  logical  unit  abs(UNIT) using format IFORM, then
 converts to multiple-precision  number  X  using  routine  MPIN.   IFORM
 should  contain  a format which allows for reading N words in A1 format,
 e.g. '(80A1)'.
0ERR  is  returned  as  .TRUE.  if  MPIN  could  not interpret input as a
 multiple-precision number or if N not positive,  otherwise  .FALSE.   In
 the former case X is returned as zero.
0Note - for a more convenient, though less flexible,  method  of  reading
 multiple-precision numbers, see MPFIN.
0MPINIT
+______
0CALL MPINIT (MP) or INITIALIZE MP
0Declares  blank  COMMON (see Section 3.1) and calls MPSET2 to initialize
 parameters.  MP is a dummy integer argument.  The Augment declaration
0      INITIALIZE MP
0causes the statement
0      CALL MPINIT (MP)
0to be generated.
0Warning  - as distributed MPINIT assumes output unit LUN = 6, at most 25
 multiple-precision digits (i.e. MT2 = 27), MNSPTR = 1, and  MXR  =  500.
 If  the  Augment  description  deck is changed, MPINIT should be changed
 accordingly (see Section 9.2).
0MPINTG
+______
0LV = MPINTG (X) or LV = INTEGR (X) or IF (INTEGR (X)) ...
0Returns .TRUE. if the (unpacked) multiple-precision number X is an exact
 integer,  .FALSE.  otherwise.  LV is a LOGICAL variable.  MPINTG must be
 declared LOGICAL if it is called directly.
0MPIO
+____
0CALL MPIO (C, N, UNIT, IFORM, ERR)
0Writes C(1), ... , C(N) in format IFORM if UNIT > 0,
 Reads  C(1), ... , C(N) in format IFORM if UNIT <= 0,
0in  both  cases  using  logical  unit  abs(UNIT).   Unformatted  I/O  is
 performed if IFORM = 'U'. C is an integer array of dimension at least N.
1Sec. 6                       MP User's Guide                     Page 48
+                             _______________
0N  and  UNIT  are  integers, ERR is a logical variable.  ERR is returned
 as .TRUE. iff N <= 0.  (It would be desirable to return ERR = .TRUE.  if
 an I/O error occurred.  This is easily done on most systems, but can not
 be done with ANSI Standard Fortran. See Section 9.3.)
0MPIS
+____
0LV = MPIS (J, K)
0Assumes that J and K are INTEGER words each containing one character  in
 the  leftmost  character  position (read under A1 format or set by a 1Hx
 data statement or unpacked by MPUPK).  MPIS returns .TRUE.  if  the  two
 characters are equal, .FALSE.  otherwise.  LV is a LOGICAL variable, and
 MPIS must be declared to be LOGICAL.
0On  most  machines,  MPIS (J, K) is equivalent to (J .EQ. K), but on the
 Burroughs B6700 it is equivalent to (J .IS. K).  For  conversion  notes,
 see Section 9.2.
0MPKADD
+______
0CALL MPKADD (X, Y, Z) or Z = X + Y
0The  same  as   CALL  MPADD  (X,  Y,  Z)  except that X and Y are packed
 multiple-precision variables.  (They may be packed or unpacked if MPKADD
 is  called  directly  -  similarly  for  MPKDIV, MPKDVI, MPKIML, MPKMLI,
 MPKMUL and MPKSUB.)  For further details, see under MPADD.
0MPKDIV
+______
0CALL MPKDIV (X, Y, Z) or Z = X/Y
0The same as  CALL MPDIV (X, Y, Z)   except  that  X  and  Y  are  packed
 multiple-precision variables.  For further details see under MPDIV.
0MPKDVI
+______
0CALL MPKDVI (X, IY, Z) or Z = X/IY
0The  same  as   CALL  MPDIVI  (X,  IY,  Z)  except  that  X  is a packed
 multiple-precision variable.  For further details, see under MPDIV.
0MPKIML
+______
0CALL MPKIML (IY, X, Z) or Z = IY*X
0The same as  CALL  MPMULI  (X,  IY,  Z)   except  that  X  is  a  packed
 multiple-precision variable.  For further details, see under MPMULI.
0MPKMLI
+______
0CALL MPKMLI (X, IY, Z) or Z = X*IY
0The  same  as   CALL  MPMULI  (X,  IY,  Z)   except  that  X is a packed
 multiple-precision variable.  For further details, see under MPMULI.
1Sec. 6                       MP User's Guide                     Page 49
+                             _______________
0MPKMUL
+______
0CALL MPKMUL (X, Y, Z) or Z = X*Y
0The same as  CALL MPMUL (X, Y,  Z)  except  that  X  and  Y  are  packed
 multiple-precision variables.  For further details, see under MPMUL.
0MPKSUB
+______
0CALL MPKSUB (X, Y, Z) or Z = X - Y
0The  same  as   CALL  MPSUB  (X,  Y,  Z)  except that X and Y are packed
 multiple-precision variables.  For further details, see under MPSUB.
0MPK3V
+_____
0CALL MPK3V (MPXX, X, Y, Z)
0Calls MPXX (XU, YU, Z) after unpacking X and Y to give XU and YU.  X and
 Y  are packed or unpacked multiple-precision variables, Z is an unpacked
 multiple-precision variable, MPXX is a subroutine taking three  unpacked
 multiple-precision  arguments  and  returning  a  result  as  the  third
 argument (e.g. MPMUL).  Called by MPKADD, MPKMUL, MPKSUB etc.
0MPK3V2
+______
0CALL MPK3V2 (MPXX, X, N, Y)
0Calls MPXX (XU, N, Z) after unpacking X to give XU.  X is  a  packed  or
 unpacked  multiple-precision variable, N is an integer, Z is an unpacked
 multiple-precision variable, MPXX is a subroutine taking three arguments
 of  the correct types and returning a result as the third argument (e.g.
 MPMULI).  Called by MPKDVI, MPKMLI etc.
0
 MPLARG
+______
0CALL MPLARG (MXINT)
0Returns  MXINT <= maximum representable INTEGER of the form  2**k  -  1.
 Integer  arithmetic  must  be performed exactly on integers in the range
 -MXINT, ... , +MXINT, so on some  machines  (e.g.  Burroughs  B6700  and
 Cyber  76) MXINT must be smaller than the largest representable integer.
 We say that the 'effective' wordlength is k+1 bits.
0Note  - The version of MPLARG supplied with the MP package does not work
 on all machines.  For conversion notes, see Section 9.2.
0MPLE
+____
0LV = MPLE (X, Y) or LV = (X .LE. Y) or IF (X .LE. Y) ...
0Returns logical value of (X < Y) for multiple-precision X and  Y.   MPLE
 must be declared type logical unless Augment interface used.  X and/or Y
 may be packed or unpacked.
1Sec. 6                       MP User's Guide                     Page 50
+                             _______________
0MPLG10
+______
0CALL MPLG10 (X, Y) or Y = LOG10 (X)
0Returns Y = ln(X)/ln(10) = log(X) to base 10, for  multiple-precision  X
 and  Y.   X  may  be  packed  or unpacked, Y is unpacked.  Uses MPLN and
 MPLNI.  For restrictions,  time  bounds,  and  the  effect  of  rounding
 options, see the description of MPLN.
0
 MPLI
+____
0CALL MPLI (X, Y) or Y = LI (X)
0Returns
0       Y = Li(X) = logarithmic integral of X
          = principal value integral from 0 to X of du/ln(u),
0using  MPEI.   X and Y are multiple-precision numbers, X > 0, X .NE.  1.
 X  may  be  packed  or  unpacked  if  MPLI is called directly.  Rounding
 options are not yet implemented.  Error in Y  could  be  induced  by  an
 O(B**(1-T))  relative perturbation in X followed by similar perturbation
 in Y.  Thus relative error in Y is small unless X is close to  1  or  to
 the zero 1.45136923488338105028...  of Li(X).  Time = O(m(T)T).
0
 MPLN
+____
0CALL MPLN (X, Y) or Y = LN (X) or Y = LOG (X) or Y = ALOG (X)
0Returns  Y  = ln(X), for multiple-precision X and Y, using MPLNS.  X may
 be packed or unpacked, Y is unpacked.  The integer part of ln(X) must be
 representable  as  a  single-precision  integer.  Time = O(sqrt(T)m(T)).
 For small integer X, MPLNI is  faster.   Asymptotically  faster  methods
 exist  (e.g.  the  Gauss-Salamin method), but are not useful unless T is
 large (see the descriptions  of  MPATAN,  MPEXP1,  MPLNGS  and  MPPIGL).
 Rounding options are implemented as for MPEXP.
0
 MPLNGM
+______
0CALL MPLNGM (X, Y) or Y = LNGM (X)
0Returns    multiple-precision    Y    =    ln(Gamma(X))   for   positive
 multiple-precision X, using Stirling's asymptotic expansion.  (X may  be
 packed  or  unpacked  if MPLNGM is called directly.)  Slower than MPGAMQ
 (unless X large) and uses more space, so use MPGAMQ and  MPLN  if  X  is
 rational and not too large, say X less than 100.  Time = O(T**3).
0Space = nl*((T+3)/2) + O(T) words, where nl is the number of terms  used
 in the asymptotic expansion, nl < 2+0.125*T*ln(B).  Thus space = O(T**2)
 words.
0Rounding options are not yet implemented, uses no guard digits.
1Sec. 6                       MP User's Guide                     Page 51
+                             _______________
0MPLNGS
+______
0CALL MPLNGS (X, Y)
0Returns   Y   =   ln(X)  for  multiple-precision  X  and  Y,  using  the
 Gauss-Salamin  algorithm  based   on   the   arithmetic-geometric   mean
 iteration.   This  is  described  in  - R. P. Brent, 'Multiple-precision
 zero-finding  methods  and  the  complexity   of   elementary   function
 evaluation',  in  Analytic Computational Complexity (ed. by J.F. Traub),
 Academic Press, 1976, 151-176.   (If abs(X-1) <  1/B,  MPLNS  is  used.)
 Time = O(ln(T)m(T)) + O(T**2) if abs(X-1) > 1/B, as for MPLNS otherwise.
0MPLNGS is slower than MPLN unless T is large (greater than  about  500),
 so  it  is  recommended for testing purposes only.  Rounding options are
 implemented as for MPEXP.
0
 MPLNI
+_____
0CALL MPLNI (N, X) or X = LNI (N) or X = ALOG (N)
                   or X = LN  (N) or X = LOG  (N)
0Returns multiple-precision X = ln(N) for small positive integer N, using
 a rapidly converging series and MPL235.  Time  =  O(T**2),  faster  than
 MPLN.  Rounding options are implemented as for MPEXP.
0
 MPLNS
+_____
0CALL MPLNS (X, Y) or Y = LNS (X)
0Returns  multiple-precision  Y  =  ln(1+X)  if X is a multiple-precision
 number satisfying the condition  abs(X) < 1/B,  error otherwise.  X  may
 be packed or unpacked if MPLNS is called directly.  Uses Newton's method
 to solve the equation exp1(-z) = X, then sets Y = -z.  (Here  exp1(z)  =
 exp(z)  - 1 is computed using MPEXP1.)  Time = O(sqrt(T)m(T)).  Rounding
 options are implemented as for MPEXP.
0
 MPLT
+____
0LV = MPLT (X, Y) or LV = (X .LT. Y) or IF (X .LT. Y) ...
0Returns logical value of (X < Y) for multiple-precision X and  Y.   MPLT
 must  be declared LOGICAL unless the Augment interface used.  X and/or Y
 may be packed or unpacked.
0
 MPL235
+______
0Computes   ln ((2**i)*(3**j)*(5**k)) = i*ln(2) + j*ln(3) + k*ln(5)   for
 small integer i, j  and  k.   ln(81/80),  ln(25/24)  and  ln(16/15)  are
 calculated  first,  using  rapidly  converging  series.  Time = O(T**2).
 Called by MPLNI, not recommended for independent use.  Rounding  options
 are not implemented, uses no guard digits.
1Sec. 6                       MP User's Guide                     Page 52
+                             _______________
0MPMAX
+_____
0CALL MPMAX (X, Y, Z) or Z = MAX (X, Y)
0Sets  Z  = max (X, Y) where X, Y and Z are multiple-precision variables.
 X and/or Y may be packed  or  unpacked  if  MPMAX  is  called  directly.
 Rounding options are irrelevant as the result is exact.
0MPMAXR
+______
0CALL MPMAXR (X) or X = CTM ('MAXR')
0Sets X to the largest possible positive multiple-precision number,  i.e.
 B**M - B**(M-T).   Rounding options are  irrelevant  as  the  result  is
 exact.
0MPMEXA
+______
0I = MPMEXA (X) or I = MAXEXP (X)
0Returns  the  maximum  allowable  exponent of multiple-precision numbers
 (the third word of COMMON /MPCOM/ - see Section  3.1).   X  is  a  dummy
 multiple-precision argument (Augment expects one).
0MPMEXB
+______
0CALL MPMEXB (I, X) or MAXEXP (X) = I
0Sets  the maximum allowable exponent of multiple-precision numbers to I.
 I should be greater than  T,  and  4*I  should  be  representable  as  a
 single-precision  integer.   X  is  a  dummy multiple-precision argument
 (Augment expects one).
0MPMIN
+_____
0CALL MPMIN (X, Y, Z) or Z = MIN (X, Y)
0Sets Z = min (X, Y) where X, Y and Z are  multiple-precision  variables.
 X  and/or  Y  may  be  packed  or  unpacked if MPMIN is called directly.
 Rounding options are irrelevant as the result is exact.
0MPMINR
+______
0CALL MPMINR (X) or X = CTM ('MINR')
0Sets X to the smallest positive  normalized  multiple-precision  number,
 i.e. X = B**(-M), where M is the  third  word  of  COMMON  /MPCOM/  (see
 Section 3.1).    Rounding options are irrelevant as the result is exact.
0MPMLP
+_____
0Performs  inner  multiplication  loop  for  MPMUL.   Carries   are   not
 propagated  in  inner  loop,  which  saves time at the expense of space.
 This routine is usually called more often than any other MP routine,  so
 it is worthwhile to optimize it.
1Sec. 6                       MP User's Guide                     Page 53
+                             _______________
0MPMOD
+_____
0CALL MPMOD (X, Y, Z) or Z = MOD (X, Y)
0Returns  Z = X - Y*int(X/Y)  for multiple-precision X, Y and Z.  X and Y
 may be both packed or both unpacked, Z is unpacked.   Here  'int'  means
 the integer part, truncated towards zero.
0Notes - 1. MPMOD returns Z as zero if Y is zero.
         2. Time = O (max (1, ln(abs(X/Y))) * T**2)),
            which is large if abs(X) >> abs(Y).
         3. Rounding options are irrelevant as the result is exact.
0
 MPMOVE
+______
0CALL MPMOVE (X, TX, Y, TY)
0Assumes  that  X  and  Y  are  multiple-precision numbers with TX and TY
 digits respectively.  X may be packed or unpacked, Y  is  unpacked  (but
 not  necessarily  initialized).   MPMOVE  moves  X  to  Y,  padding with
 trailing zeros if TX < TY, or rounding as specified by  RNDRL  (see  the
 description  of  MPCQM)  if  TX  >  TY.  Note that TX and TY are integer
 arguments, and the value  of  T  (second  word  of  COMMON  /MPCOM/)  is
 irrelevant.
0MPSTR (X, Y) is equivalent and faster if TX = TY = T and X is  unpacked.
0
0MPMUL
+_____
0CALL MPMUL (X, Y, Z) or Z = X*Y
0Multiplies  X  and  Y,  returning  result  in  Z, for multiple-precision
 variables X, Y and Z.  The simple O(T**2) algorithm  is  used,  with  at
 least two guard digits.  Advantage is taken of any zero digits in X, but
 not in Y.  Asymptotically faster algorithms are known (see  for  example
 Knuth,  Vol.  2),  but  are  difficult  to  implement  in  Fortran in an
 efficient and machine-independent manner.
0In description of other multiple-precision routines, m(T) is the time to
 perform T-digit multiple-precision multiplication.  Thus m(T) =  O(T**2)
 with  the present version of MPMUL, but m(T) = O(T.log(T)log(log(T))) is
 theoretically possible (see Knuth, Vol. 2).
0Rounding  options are implemented.  They are controlled by the parameter
 RNDRL in COMMON /MPCOM/ (see Section 3.1) as follows -
0  RNDRL = 0 - truncate towards zero, error less than 1.01 units  in  the
               last place (ulp),
0  RNDRL = 1 - round to nearest representable multiple-precision  number,
               preferring last digit even in case of a tie, error at most
               0.5 ulp (unless the result underflows),
1Sec. 6                       MP User's Guide                     Page 54
+                             _______________
0  RNDRL = 2 - round down (towards -infinity), error less than 1 ulp,
0  RNDRL = 3 - round up (towards +infinity), error less than 1 ulp.
0Note that RNDRL = 0 is the fastest (approximately twice as fast  as  the
 other  options) because less guard digits are required.  Thus, RNDRL = 0
 is recommended for general use and is the default option.
0
 MPMULI
+______
0CALL MPMULI (X, IY, Z) or Z = X*IY
0Multiplies multiple-precision X by single-precision  integer  IY  giving
 multiple-precision  Z.   Time  =  O(T),  which  is  faster  than  MPMUL.
 Multiplication by 1 may be used to normalize a number even if  the  last
 digit is B.  Rounding options are implemented as for MPCQM.
0
 MPMULQ
+______
0CALL MPMULQ (X, I, J, Y) or Y = MULQ (X, I, J)
0Multiplies   multiple-precision   X   by  rational  number  I/J,  giving
 multiple-precision Y.  Here I and J  are  single-precision  integers,  J
 nonzero.  Time = O(T).
0Rounding options are determined by RNDRL in COMMON /MPCOM/ (see  Section
 3.1) as follows -
0  RNDRL = 0 - as for MPMULI followed by MPDIVI,
0  RNDRL = 1, 2, or 3 - as for MPCQM.
0Augment users should note that Y = MULQ (X, I, J) is usually faster than
 Y = X * CTM (I,J).
0
 MPMULS
+______
0CALL MPMULS (X, I, J, K, L)
0Sets  X = X*I*J/(K*L)  for unpacked multiple-precision X, integer I,  J,
 K  and  L  (K*L  nonzero).  Calls MPMULQ once if I*J and K*L are not too
 large, otherwise calls MPMULQ twice.  Rounding is not best possible, but
 the directed rounding options give correct upper and lower bounds.
0
 MPNE
+____
0LV = MPNE (X, Y) or LV = (X .NE. Y) or IF (X .NE. Y) ...
0Returns  logical  value  of  (X  .NE. Y) for multiple-precision X and Y.
 MPNE must be declared LOGICAL unless  the  Augment  interface  used.   X
 and/or Y may be packed or unpacked.
1Sec. 6                       MP User's Guide                     Page 55
+                             _______________
0MPNEG
+_____
0CALL MPNEG (X, Y) or Y = -X
0Sets Y = -X for multiple-precision numbers X and Y.  X and Y may be both
 packed or both unpacked.  (This violates the usual rule of returning  an
 unpacked  result  regardless  of  whether  the  argument  is  packed  or
 unpacked.  The reason is that Augment does not allow unary operators  to
 return  a  type  different  from that of the variable which they operate
 on.)
0
0MPNEW
+_____
0CALL MPNEW (I)
0Returns  integer  index I such that words I, ... , I+T+1 of blank COMMON
 are available for use.  Sets  SPTR  =  I+T+2  and  updates  MXSPTR  (see
 Section 3.1).  MPSTOV is called if MXR is too small.
0Note - CALL MPNEW (I) is equivalent to CALL MPNEW2 (I, T+2).
0
 MPNEW2
+______
0CALL MPNEW2 (I, J)
0Returns integer index I such that words I, ...  ,  I+abs(J)-1  of  blank
 COMMON  are  available for use.  Sets SPTR = I+abs(J) and updates MXSPTR
 (see Section 3.1).  MPSTOV is called if MXR is too small.
0
 MPNZR
+_____
0Normalizes  and  rounds  a  multiple-precision  number (represented in a
 nonstandard format, with separate sign and exponent, and  possibly  some
 guard  digits).   Rounding  depends on RNDRL in COMMON /MPCOM/ - see the
 description of MPCQM.  Called by MPADD, MPDIVI, MPMUL, MPMULI, etc., and
 not recommended for independent use.
0
 MPOUT
+_____
0CALL MPOUT (X, C, P, N)
0Converts  multiple-precision X to FP.N format in C, which may be printed
 under PA1 format.  (Here FP.N means sign, P-N-2 'decimal' places  before
 the  'decimal'  point,  and  N  places  after  it.)  Note that N = -1 is
 allowed, and effectively gives IP format (i.e., sign and  P-1  'decimal'
 digits).  Digits after the 'decimal' point are blanked out if they could
 not be significant.  The output base is OUTBAS in  COMMON  /MPCOM/  (see
 Section  3.1),  in  the range two to sixteen, with default value ten.  P
 and N are integers, C is an integer array of dimension at  least  P.   X
 may be packed or unpacked.
1Sec. 6                       MP User's Guide                     Page 56
+                             _______________
0Rounding is defined  by  the  parameter  RNDRL  in  COMMON  /MPCOM/,  as
 follows -
0  0 or 1 - Round   to   approximately   the   nearest   decimal   number
            representable in the specified  format.   The  error  in  the
            conversion  to  decimal  is  less  than 0.6 units in the last
            decimal place.
0  2 -      Round   down   (towards   -infinity)   to  a  decimal  number
            representable in the specified format.
0  3 -      Round   up   (towards   +infinity)   to   a   decimal  number
            representable in the specified format.
0Efficiency is higher if B is a power of 10 than if not.  Time = O(T**2).
0Note  -  MPFOUT  and  MP40D  are  more  convenient to use, although less
 flexible than MPOUT.
0
 MPOUTE
+______
0CALL MPOUTE (X, C, J, P)
0Assumes X is a multiple-precision number  and  C  an  integer  array  of
 dimension  at  least  P  > 3.  J and P are integers.  On return J is the
 exponent (to base OUTBAS) of X and the fraction is in  C,  ready  to  be
 printed  in  A1  format.   For example, we could print J and C in format
 '(I10, 1X, PA1)'.  The fraction has one place before the 'decimal' point
 and P-3 after.  The output base is OUTBAS in COMMON /MPCOM/ (see Section
 3.1), with default value ten.  Rounding is not the  best  possible,  but
 the  directed  rounding options (RNDRL = 2 and 3) give correct lower and
 upper bounds on the true result (see the description of MPOUT).   X  may
 be packed or unpacked.
0Note - MPFOUT  is  a  more  convenient  (though  less  flexible)  output
 routine.
0
 MPOUTF
+______
0CALL MPOUTF (X, P, N, IFORM, ERR)  or  ERR = MPOUTF (X, P, N, IFORM)
                                    or  IF (MPOUTF (X, P, N, IFORM))  ...
0Writes multiple-precision number X on logical unit LUN (fourth  word  of
 COMMON  /MPCOM/)  in  format  IFORM  after  converting to FP.N 'decimal'
 representation using MPOUT.  The output base  is  OUTBAS,  with  default
 value  ten.  For further details see description of MPOUT.  IFORM should
 contain format which allows for output of P words in A1 format, plus any
 desired  spacing  etc.   ERR  is  returned  as  .TRUE.  iff  P <= 0 (see
 description of MPIO).  X may be packed or unpacked if MPFOUT  is  called
 directly.  For rounding options see MPOUT.  Time = O(T**2).
0Note - a more  convenient  (though  less  flexible)  output  routine  is
 MPFOUT.
1Sec. 6                       MP User's Guide                     Page 57
+                             _______________
0MPOUT2
+______
0CALL MPOUT2 (X, C, P, N, NB)
0Same as MPOUT except that output representation is in base NB, where two
 < NB < sixteen, e.g. NB = eight gives octal output, NB =  sixteen  gives
 hexadecimal.     Output   digits   are   0123456789ABCDEF.    X   is   a
 multiple-precision number, P, N and NB are integers,  C  is  an  integer
 array  of  dimension  at least P.  Time = O(T**2).  For rounding options
 see the description of MPOUT.
0Note  - MPOUT2 is superfluous now that OUTBAS determines the output base
 for MPOUT.  It is included for compatibility with  earlier  versions  of
 the MP package.
0
 MPOVFL
+______
0CALL MPOVFL (X)
0Called  on multiple-precision overflow, that is when the exponent of the
 multiple-precision number X would exceed  M.   Execution  is  terminated
 with  an  error message after calling MPMAXR(X) and setting MXEXPN = M+1
 (see Section 3.1).
0
 MPPACK
+______
0CALL MPPACK (X, Y) or Y = X or Y = CTP (X)
0Assumes that X is a multiple-precision number represented as usual in an
 integer  array  of  dimension at least T+2, and Y is an integer array of
 dimension at least int((T+3)/2).  X is stored in a compact format in  Y,
 and may be retrieved by calling MPUNPK (Y, X).
0MPPACK and MPUNPK are useful if space  is  critical,  for  example  when
 working  with  large  arrays  of  multiple-precision  numbers.  X may be
 packed or unpacked if MPPACK is called directly.  (If X  is  packed  the
 effect of MPPACK is the same as that of MPSTR.)
0Augment users - X is type MULTIPLE, Y is type MULTIPAK.
0
 The packed format is as follows -
0      word 2 = exponent (as for unpacked multiple-precision numbers),
0      words 1 and 3, ... , int((T+3)/2) = base B digits packed  two  per
                word  (digits  pq  are  represented  as p*B+q), with sign
                attached to word 1.
0Thus, zero is represented in both packed and unpacked formats by a  zero
 first  word,  with  the  following  words undefined.  The first word has
 absolute value greater than 1 if  and  only  if  the  multiple-precision
 number is nonzero and is represented in the packed format.
1Sec. 6                       MP User's Guide                     Page 58
+                             _______________
0MPPARA
+______
0I = MPPARA (A) or I = PARAM (A)
0Returns the integer value of the word of COMMON /MPCOM/ corresponding to
 the Hollerith string  A.   For  the  allowable  values  of  A,  see  the
 description of MPPARM.
0
 MPPARB
+______
0CALL MPPARB (I, A) or PARAM (A) = I
0Sets  the word of COMMON /MPCOM/ corresponding to the Hollerith string A
 to the integer value  I.   For  the  allowable  values  of  A,  see  the
 description  of MPPARM.  For the allowable values of I, see Section 3.1.
0
 MPPARC
+______
0CALL MPPARC (N, J)
0Sets the N-th word of COMMON /MPCOM/ to the integer value J, for
 0 < N < 24.  Faster than MPPARB, but less mnemonic.
0
0MPPARM
+______
0CALL MPPARM (HOLL, SET, J)
0HOLL is a Hollerith string, SET is a  Boolean  variable,  and  J  is  an
 integer.  The first three characters of HOLL should agree with the first
 three characters of one of the following -
0  (1)  'BASE'        (2)  'NUMDIG'      (3)  'MAXEXP'      (4)  'LUN'
   (5)  'MXR'         (6)  'SPTR'        (7)  'MXSPTR'      (8)  'MNSPTR'
   (9)  'MXEXPN'      (10) 'MNEXPN'      (11) 'RNDRL'       (12) 'KTUNFL'
   (13) 'MXUNFL'      (14) 'DECPL'       (15) 'MT2'         (16) 'MXINT'
   (17) 'EXWID'       (18) 'INRECL'      (19) 'INBASE'      (20) 'OUTBAS'
   (21) 'EXPCH'       (22) 'CHWORD'      (23) 'ONESCP'
0The corresponding word in COMMON /MPCOM/ will be set to the integer J if
 SET = .TRUE., or its value will be returned in J if SET =  .FALSE.   For
 the meanings and default settings of these parameters in COMMON /MPCOM/,
 see Section 3.1.
0
0MPPARN
+______
0J = MPPARN (N)
0Returns  the  value  of the N-th word of COMMON /MPCOM/, for 0 < N < 24,
 zero otherwise.  Faster than MPPARA, but less mnemonic.
1Sec. 6                       MP User's Guide                     Page 59
+                             _______________
0MPPI
+____
0
 CALL MPPI (X) or X = CTM ('PI')
0Sets multiple-precision X = pi = 3.14159... to the available  precision.
 Uses the formula
0      pi = 16*arctan(1/5) - 4*arctan(1/239)
0of Machin.  Time = O(T**2).  (For asymptotically faster methods, see the
 description of MPPIGL.)  Rounding options are implemented as for  MPEXP.
0
0MPPIGL
+______
0CALL MPPIGL (X)
0Sets  multiple-precision  X  =  pi  = 3.14159... to almost the available
 precision.  Uses the Gauss-Legendre  algorithm.   This  method  requires
 time  O(ln(T)m(T)),  so  it  is  slower than MPPI if m(T) = O(T**2), but
 would be faster for large T if a faster  multiplication  algorithm  were
 used  (see  the  description  of  MPMUL).   For  a  description  of  the
 Gauss-Legendre  algorithm  see  -  R.  P.   Brent,   'Multiple-precision
 zero-finding  and  the complexity of elementary function evaluation', in
 Analytic Computational  Complexity  (edited  by  J.F.  Traub),  Academic
 Press, l976, 151-176.
0Rounding options are not implemented, uses no guard digits.  Recommended
 for testing purposes only.
0
0MPPOLY
+______
0
 CALL MPPOLY (X, Y, IC, N)
0Sets Y = IC(1) + IC(2)*X + ... +  IC(N)*X**(N-1),  where  X  and  Y  are
 multiple-precision  numbers,  and IC is an integer array of dimension at
 least N > 0.  Rounding is not best possible, but the  directed  rounding
 options  (RNDRL  =  2  and 3) give correct lower and upper bounds on the
 true result.
0
0MPPWR
+_____
0
 CALL MPPWR (X, N, Y) or Y = X**N
0Returns  Y = X**N,  for  multiple-precision  X and Y,  integer  N,  with
 0**0 = 1.  X may be packed or unpacked, Y is unpacked.  Rounding options
 are implemented as for MPEXP.
1Sec. 6                       MP User's Guide                     Page 60
+                             _______________
0MPPWRA
+______
0CALL MPPWRA (X, N, Y)
0Returns  Y =  X**abs(N)  for  (unpacked)  multiple-precision  X  and  Y,
 integer   N,  with  0**0  =  1.   Uses  truncated  rather  than  rounded
 arithmetic, and no guard digits,  but  directed  rounding  options  give
 correct  upper  and  lower  bounds.   Called  by  MPEXP, MPPWR etc., not
 recommended for independent use (use MPPWR instead).
0
 MPPWR2
+______
0CALL MPPWR2 (X, Y, Z) or Z = X**Y
0Returns Z = X**Y for multiple-precision numbers X, Y and Z, where  X  is
 positive  (X  =  0  is allowed if Y > 0).  X and Y may be both packed or
 both unpacked, Z is unpacked.  Uses MPLN and MPEXP, so slower than MPPWR
 and MPQPWR.  Rounding options are implemented as for MPEXP.
0
 MPQPWR
+______
0CALL MPQPWR (I, J, K, L, X) or X = QPWR (I, J, K, L)
0Sets  multiple-precision  X  =  (I/J)**(K/L)  for integers I, J, K and L
 (with the obvious conditions to ensure that a real result  is  defined).
 Uses  MPROOT  if  abs(L)  is  small,  otherwise  uses  MPLNI  and MPEXP.
 Rounding options are implemented as for MPEXP.
0Augment users - X = QPWR (I, J, K, L) is usually faster than
                 X = CTM (I, J) ** CTM (K, L).
0
0MPREC
+_____
0CALL MPREC (X, Y) or Y = REC (X)
0Returns  Y  =  1/X, for multiple-precision X and Y.  Uses MPDIVL if T is
 small or RNDRL > 0, MPROOT otherwise.  Time =  O(m(T))  if  RNDRL  =  0,
 O(T**2) if RNDRL > 0.
0Augment users - Y = REC (X) is faster than Y = 1/X.
0
0MPRESN
+______
0CALL MPRESN (N)
0Restores  T,  M  and  RNDRL  which are assumed to have been saved on the
 stack in blank COMMON by a call to MPSAVN (N).  The integer N should  be
 as returned by MPSAVN.  SPTR is restor