! Copyright (c) 1993 Unicomp, Inc.
! The modifications made by Unicomp, Inc. are not to be
! distributed or copied without permission from Unicomp, Inc.

MODULE XP_REAL_MODULE

!  The programs in this module are based on the following work:

!      ALGORITHM 693, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 17, NO. 2, PP. 273-283.  JUNE, 1991.

!  The Fortran 77 version was written by David M. Smith as:

!     FM 1.0                  David M Smith                  1-06-90

!  The Fortran 90 interface was written by Walt Brainerd, Unicomp, Inc.

!  The XP_REAL package performs floating point extended precision arithmetic.


!  Subroutine XP_SET sets the working precision, which is the value of the
!  variable NPREC; its argument is the number of decimal digits of precision.

!  JBASE is the base in which the arithmetic is done.
!        JBASE must be bigger than one, and less than or equal
!        to the square root of the largest representable integer.
!        For best efficiency JBASE should be about 1/4 to 1/2 as big as
!        the square root of the largest representable integer.
!        Input and output conversion are much faster when JBASE is a
!        power of ten.

!  NDIG  is the number of base JBASE digits which are carried in the
!        multiple precision numbers.  NDIG must be at least two.
!        The upper limit for NDIG is defined in the PARAMETER statement
!        at the top of each routine and is restricted only by the amount
!        of memory available.

!  The format of a floating multiple precision number is as
!  follows.  The number is kept in an integer array with the first
!  element containing the exponent and each of the succeeding NDIG
!  locations containing one digit of the mantissa, expressed in base
!  JBASE.  The exponent is a power of JBASE and the implied radix point
!  is immediately before the first digit of the mantissa.  Every nonzero
!  number is normalized so that the second array element (the first
!  digit of the mantissa) is nonzero.

!  The sign of the number is carried on the second array element only.
!  Elements 3,4,... are always nonnegative.
!  The exponent is a signed integer and may be as large in magnitude as
!  MXEXP (defined in XP_SET).

!  For JBASE = 10,000 and NDIG = 4 the number -pi would have this
!  representation:
!                   Word 1         2         3         4         5

!                        1    -31415  92653590

!  Because of normalization, the equivalent number of base 10
!  significant digits for an XP number may be as small as
!  LOG10(JBASE)*(NDIG-1) + 1.


!  Subroutine FMOUT performs conversion of an FM number to base 10 and
!  formats it for output as a character array.
!  The user sets JFORM1 and JFORM2 to determine the output format.

!  JFORM1 = 0     E  format       ( .314159M+6 )
!         = 1     ES format       ( 3.14159M+5 )
!         = 2     F  format       ( 314159.000 )

!  JFORM2 = Number of significant digits to display (if JFORM1 = 0, 1).
!           If JFORM2.EQ.0 then a default number of digits is chosen.
!           The default is roughly the full precision of the number.
!  JFORM2 = Number of digits after the decimal point (if JFORM1 = 2).
!           See the FMOUT documentation for more details.


!  KW is the unit number to be used for all output from the package,
!     including error and warning messages, and trace output.


!  NTRACE and LVLTRC control trace printout from the package.

!  NTRACE =  0   No printout except warnings and errors.
!         =  1   The result of each call to one of the routines
!                   is printed in base 10, using FMOUT.
!         = -1   The result of each call to one of the routines
!                   is printed in internal base JBASE format.
!         =  2   The input arguments and result of each call to one
!                   of the routines is printed in base 10, using FMOUT.
!         = -2   The input arguments and result of each call to one
!                   of the routines is printed in base JBASE format.


!  LVLTRC defines the call level to which the trace is done.  LVLTRC = 1
!         means only FM routines called directly by the user are traced,
!         LVLTRC = 2 also prints traces for FM routines called by other
!         FM routines called directly by the user, etc.

!  In the above description internal JBASE format means the number is
!  printed as it appears in the array as an exponent followed by NDIG
!  base JBASE digits.

!  KFLAG is a condition parameter returned by the package.  Negative
!        values indicate conditions for which a warning message will
!        be printed unless KWARN = 0.  Positive values indicate
!        conditions which may be of interest but are not errors.
!        No warning message is printed if KFLAG is nonnegative.

!    KFLAG =  0     Normal operation.

!          =  1     One of the operands in FMADD or FMSUB was
!                       insignificant with respect to the other, so
!                       that the result was equal to the argument of
!                       larger magnitude.
!          =  2     In converting an FM number to a one word integer
!                       in FMM2I the FM number was not exactly an
!                       integer.  The next integer toward zero was
!                       returned.

!          = -1     NDIG was less than 2 or more than MXNDIG.
!          = -2     JBASE was less than 2 or more than MXBASE.
!          = -3     An exponent was out of range.
!          = -4     Invalid input argument(s) to an FM routine.
!                        UNKNOWN was returned.
!          = -5     + or - OVERFLOW was generated as a result from an
!                        FM routine.
!          = -6     + or - UNDERFLOW was generated as a result from an
!                        FM routine.
!          = -7     The input string to FMINP was not legal.
!          = -8     The output array for FMOUT was not large enough.
!          = -9     Precision could not be raised enough to provide all
!                        requested guard digits.  UNKNOWN was returned.
!          = -10    An FM input argument was too small in magnitude to
!                        convert in FMM2SP or FMM2DP.  Zero was
!                        returned.

!  When a negative KFLAG condition is encountered the routine calls
!  FMWARN, which uses the value of KWARN as follows.

!  KWARN = 0     Execution continues and no message is printed.
!        = 1     A warning message is printed and execution continues.
!        = 2     A warning message is printed and execution stops.

!  When an overflow or underflow is generated for an operation in which
!  an input argument was already an overflow or underflow, no additional
!  message is printed.  When an unknown result is generated and an input
!  argument was already unknown, no additional message is printed.  In
!  these cases the negative KFLAG value is still returned.


!  KRAD = 0     Causes all angles in the trigonometric functions and
!                  inverse functions to be measured in degrees.
!       = 1     Causes all angles to be measured in radians.


!  KROUND = 0   Causes all final results to be chopped (rounded toward
!                  zero).  Intermediate results are rounded.
!         = 1   Causes all results to be rounded to the nearest FM
!                  number, or to the value with an even last digit if
!                  the result is halfway between two FM numbers.


!  Here is a list of the routines in XP that are designed to be called
!  by the user.

   INTERFACE ASSIGNMENT (=)
      MODULE PROCEDURE XP_GETS_INT, XP_GETS_REAL, XP_GETS_DP
   END INTERFACE

   INTERFACE ABS
      MODULE PROCEDURE XP_ABS
   END INTERFACE

   INTERFACE ACOS
      MODULE PROCEDURE XP_ACOS
   END INTERFACE

   INTERFACE OPERATOR (+)
      MODULE PROCEDURE XP_ADD, XPI_ADD, IXP_ADD
   END INTERFACE

   INTERFACE AINT
      MODULE PROCEDURE XP_AINT
   END INTERFACE

   INTERFACE ASIN
      MODULE PROCEDURE XP_ASIN
   END INTERFACE

   INTERFACE ATAN
      MODULE PROCEDURE XP_ATAN
   END INTERFACE

   INTERFACE ATAN2
      MODULE PROCEDURE XP_ATAN2
   END INTERFACE

   INTERFACE HUGE
      MODULE PROCEDURE XP_HUGE
   END INTERFACE

   INTERFACE OPERATOR (==)
      MODULE PROCEDURE XP_EQ
   END INTERFACE

   INTERFACE OPERATOR (/=)
      MODULE PROCEDURE XP_NE
   END INTERFACE

   INTERFACE OPERATOR (<=)
      MODULE PROCEDURE XP_LE
   END INTERFACE

   INTERFACE OPERATOR (>=)
      MODULE PROCEDURE XP_GE
   END INTERFACE

   INTERFACE OPERATOR (<)
      MODULE PROCEDURE XP_LT
   END INTERFACE

   INTERFACE OPERATOR (>)
      MODULE PROCEDURE XP_GT
   END INTERFACE

   INTERFACE COS
      MODULE PROCEDURE XP_COS
   END INTERFACE

   INTERFACE COSH
      MODULE PROCEDURE XP_COSH
   END INTERFACE

   INTERFACE DIM
      MODULE PROCEDURE XP_DIM
   END INTERFACE

   INTERFACE OPERATOR (/)
      MODULE PROCEDURE IXP_DIV, XPI_DIV, XP_DIV
   END INTERFACE

   INTERFACE XP
      MODULE PROCEDURE XP_I2M, XP_R2M, XP_DP2M, XP_C2M
   END INTERFACE

   INTERFACE EXP
      MODULE PROCEDURE XP_EXP
   END INTERFACE

   INTERFACE OPERATOR (**)
      MODULE PROCEDURE IXP_PWR, XPI_PWR, XP_PWR
   END INTERFACE

   INTERFACE LOG10
      MODULE PROCEDURE XP_LOG10
   END INTERFACE

   INTERFACE LOG
      MODULE PROCEDURE XP_LOG
   END INTERFACE

   INTERFACE DBLE
      MODULE PROCEDURE XP_M2DP
   END INTERFACE

   INTERFACE INT
      MODULE PROCEDURE XP_M2I
   END INTERFACE

   INTERFACE REAL
      MODULE PROCEDURE XP_M2SP
   END INTERFACE

   INTERFACE MAX
      MODULE PROCEDURE XP_MAX
   END INTERFACE

   INTERFACE MIN
      MODULE PROCEDURE XP_MIN
   END INTERFACE

   INTERFACE MOD
      MODULE PROCEDURE XP_MOD
   END INTERFACE

   INTERFACE OPERATOR (*)
      MODULE PROCEDURE XP_MPY, XPI_MPY, IXP_MPY
   END INTERFACE

   INTERFACE ANINT
      MODULE PROCEDURE XP_ANINT
   END INTERFACE

   INTERFACE PI
      MODULE PROCEDURE XP_PI
   END INTERFACE

   INTERFACE SIGN
      MODULE PROCEDURE XP_SIGN
   END INTERFACE

   INTERFACE SIN
      MODULE PROCEDURE XP_SIN
   END INTERFACE

   INTERFACE SINH
      MODULE PROCEDURE XP_SINH
   END INTERFACE

   INTERFACE SQRT
      MODULE PROCEDURE XP_SQRT
   END INTERFACE

   INTERFACE OPERATOR (-)
      MODULE PROCEDURE XP_SUB
   END INTERFACE

   INTERFACE TAN
      MODULE PROCEDURE XP_TAN
   END INTERFACE

   INTERFACE TANH
      MODULE PROCEDURE XP_TANH
   END INTERFACE

!   INTERFACE PRINT
!      MODULE PROCEDURE XP_PRINT
!   END INTERFACE

!   PRIVATE XP_PRINT


!  FM numbers in unpacked format have size LUNPCK.


!  PORTABILITY NOTES.

!  In routines FMEQU and FMSUB there is code which attempts to
!  determine if two input arrays refer to identical memory locations.
!  Some optimizing compilers assume the arrays must be distinct and
!  may remove code which would then be redundant.  This code removal
!  could cause errors, so the tests are done in a way which should
!  keep the compiler from removing code. The current version works
!  correctly on all compilers tested.  Computing SIN(1.0) in radian
!  mode should reveal whether other compilers handle it correctly.
!  If there is a problem, SIN(1) gives 0.999... instead of 0.841....
!  To fix such a problem, MB can be copied to a local temporary array
!  and then negated in FMSUB before calling FMADD2.  For FMEQU
!  simply set KSAME = 0 after the section which tries to determine if
!  MA and MB are the same array. This makes both routines run slower.
!  A simpler fix which often works is to re-compile at a lower
!  optimization (e.g., OPT=0).

!  Using the CFT77 compiler on a Cray X-MP computer there are
!  problems using a large value for JBASE when the default 46-bit
!  integer arithmetic mode is used.  In particular, FMSET chooses
!  too large a JBASE value since some of the arithmetic in the
!  MAXINT loop is done with 64-bit precision.  Setting JBASE = 10**6
!  or less may be ok, but the preferred solution is to select the
!  64-bit integer compiler option.  Then JBASE = 10**9 can be used.

! --------------------------------------------------------------------

!  The size of all arrays is controlled by defining two parameters:
!  MXNDIG is the maximum value the user can set NDIG,
!  NBITS  is the number of bits per integer word.

      PARAMETER ( MXNDIG=256 , NBITS=32 ,  &

!             Define the array sizes:

                LUNPCK = (6*MXNDIG)/5 + 20,  &
                LMWA   = 2*LUNPCK         , LJSUMS = 8*LUNPCK,  &
                LMBUFF = ((LUNPCK+3)*(NBITS-1)*301)/2000 + 6  )

   TYPE XP_REAL
      INTEGER, DIMENSION (LUNPCK) :: DIGITS
   END TYPE XP_REAL

!  The following are work arrays used by the low-level arithmetic routines,
!  definitions for overflow and underflow thresholds, etc.


      INTEGER, PRIVATE :: MWA(LMWA), KSTACK(19)

!  The following contain information about saved constants.

      INTEGER, PRIVATE ::   &
                          MPISAV(LUNPCK),MESAV(LUNPCK),MLBSAV(LUNPCK),  &
                          MLN1(LUNPCK),MLN2(LUNPCK),MLN3(LUNPCK),  &
                          MLN4(LUNPCK)

!             MJSUMS is an array which can contain several XP numbers
!             being used to accumulate the concurrent sums in FMEXP2
!             and FMSIN2.  When MXNDIG is 256, eight is about the maximum
!             number of sums needed (but this depends on JBASE).  For
!             larger MXNDIG dimensioning MJSUMS to hold more than eight
!             FM numbers could speed the functions up.

      INTEGER, PRIVATE :: MJSUMS(LJSUMS)

!             MBUFF is a character array used by FMPRNT for printing
!             output from FMOUT.  This array may also be used for calls
!             to FMOUT from outside the XP package.

      CHARACTER, PRIVATE :: MBUFF(LMBUFF)

!             The following are scratch arrays for temporary storage of XP
!             numbers while computing various functions.

      INTEGER, PRIVATE :: M01(LUNPCK),M02(LUNPCK),M03(LUNPCK),M04(LUNPCK),  &
                          M05(LUNPCK),M06(LUNPCK)

!             The following are scratch arrays used to hold input arguments
!             in unpacked format.

!  Base and precision will be set to give at least NPREC+3 decimal
!  digits of precision (giving the user three base ten guard digits).

!  JBASE is set to the largest permissible power of ten.
!  JFORM1 and JFORM2 are set to ES format displaying NPREC significant
!  digits.

!  The trace option is set off, and the mode for angles in trig
!  functions is set to radians.  The rounding mode is set to symmetric
!  rounding.

!  KW, the unit number for all XP output, is set to 6.
!  This includes trace output and error messages.

   INTEGER :: KW = 6

   INTEGER, PRIVATE, PARAMETER :: MAXINT = HUGE (0)

   DOUBLE PRECISION, PRIVATE, PARAMETER :: DPMAX = HUGE (0.0D0)

   REAL, PRIVATE, PARAMETER :: SPMAX = HUGE (0.0) / 1.01

!             MXNDG2 is the maximum value for NDIG which can be used
!                    internally.  Several of the FM routines may raise
!                    NDIG above MXNDIG temporarily, in order to
!                    compute correctly rounded results.
!                    In the definition of LUNPCK The '6/5' condition
!                    allows for converting from a large base to the
!                    (smaller) largest power of ten base for output
!                    conversion.
!                    The '+ 20' condition allows for the need to carry
!                    many guard digits when using a small base like 2.

   INTEGER, PRIVATE, PARAMETER :: MXNDG2 = LUNPCK - 1

!             MXBASE is the maximum value for JBASE.
!             JBASE is the currently used base for arithmetic.

   INTEGER, PRIVATE, PARAMETER :: K_RANGE = RANGE (0) / 2
   INTEGER, PRIVATE, PARAMETER :: MXBASE = 10 ** K_RANGE
   INTEGER :: JBASE = MXBASE

!             NDIG is the number of digits currently being carried.
!             It is initially set to 40.

   INTEGER, PARAMETER, PRIVATE :: INITIAL_NPREC = 40
   INTEGER :: NPREC = INITIAL_NPREC
   INTEGER, PRIVATE :: NPSAVE = INITIAL_NPREC
!   INTEGER, PRIVATE :: NDIG = 2 + (INITIAL_NPREC+2)/K_RANGE
   INTEGER, PUBLIC :: NDIG = 2 + (INITIAL_NPREC+2)/K_RANGE

!             KFLAG is the flag for error conditions.

   INTEGER :: KFLAG = 0

!             NTRACE is the trace switch.  Default is no printing.

   INTEGER :: NTRACE = 0

!             LVLTRC is the trace level.  Default is to trace only
!             routines called directly by the user.

  INTEGER ::  LVLTRC = 1

!             NCALL is the call stack pointer.

  INTEGER, PRIVATE ::  NCALL = 0

!             Some constants which are often needed are stored in the
!             maximum precision which they have been computed with the
!             currently used base.  This speeds up the trig, log, power,
!             and exponential functions.

!             NDIGPI is the number of digits available in the currently
!                    stored value of pi (MPISAV).

  INTEGER, PRIVATE :: NDIGPI = 0

!             NJBPI is the value of JBASE for the currently stored
!                   value of pi.

  INTEGER, PRIVATE :: NJBPI = 0

!             NDIGE is the number of digits available in the currently
!                   stored value of e (MESAV).

  INTEGER, PRIVATE :: NDIGE = 0

!             NJBE is the value of JBASE for the currently stored
!                  value of e.

  INTEGER, PRIVATE :: NJBE = 0

!             NDIGLB is the number of digits available in the currently
!                     stored value of LN(JBASE) (MLBSAV).

  INTEGER, PRIVATE :: NDIGLB = 0

!             NJBLB is the value of JBASE for the currently stored
!                   value of LN(JBASE).

  INTEGER, PRIVATE :: NJBLB = 0

!             NDIGLI is the number of digits available in the currently
!                    stored values of the four logarithms
!                    used by FMLNI: MLN1 - MLN4.

  INTEGER, PRIVATE :: NDIGLI = 0

!             NJBLI is the value of JBASE for the currently stored
!                   values of MLN1 - MLN4.

  INTEGER, PRIVATE :: NJBLI = 0

!             MXEXP  is the current maximum exponent.
!             MXEXP2 is the internal maximum exponent. This is used to
!                    define the overflow and underflow thresholds.

!             These values are chosen so that FM routines can raise the
!             overflow/underflow limit temporarily while computing
!             intermediate results, and so that EXP(MAXINT) is greater
!             than MXBASE**(MXEXP2+1).

!             The overflow threshold is JBASE**(MXEXP+1), and the
!             underflow threshold is JBASE**(-MXEXP-1).
!             This means the valid exponents in the first word of an FM
!             number can range from -MXEXP to MXEXP+1 (inclusive).

  INTEGER, PRIVATE, PARAMETER :: MXEXP_PARAM = (MXBASE*MXBASE)/(4.6051702*K_RANGE) - 1
  INTEGER, PRIVATE :: MXEXP = MXEXP_PARAM
  INTEGER, PRIVATE, PARAMETER :: MXEXP2 = 2*MXEXP_PARAM + MXEXP_PARAM/100

!             KARGSW is a switch used by some of the elementary function
!                    routines to disable argument checking while doing
!                    calculations where no exceptions can occur.
!                    See FMARGS for a description of argument checking.
!                    KARGSW = 0 is the normal setting,
!                           = 1 means argument checking is disabled.

  INTEGER, PRIVATE :: KARGSW = 0

!             KEXPUN is the exponent used as a special symbol for
!                    underflowed results.

  INTEGER, PRIVATE, PARAMETER :: KEXPUN = -MXEXP2 - 5*MXNDIG

!             KEXPOV is the exponent used as a special symbol for
!                    overflowed results.

  INTEGER, PRIVATE, PARAMETER :: KEXPOV = -KEXPUN

!             KUNKNO is the exponent used as a special symbol for
!                    unknown FM results (1/0, SQRT(-3.0), etc).

  INTEGER, PRIVATE, PARAMETER :: KUNKNO = KEXPOV + 5*MXNDIG

!             RUNKNO is returned from FM to real or double conversion
!                    routines when no valid result can be expressed in
!                    real or double precision.  On systems which provide
!                    a value for undefined results (e.g., Not A Number)
!                    setting RUNKNO to that value is reasonable.  On
!                    other systems set it to a value which is likely to
!                    make any subsequent results obviously wrong which
!                    use it.  In either case a KFLAG = -4 condition is
!                    also returned.

   REAL, PRIVATE, PARAMETER :: RUNKNO = -1.01*SPMAX

!             IUNKNO is returned from FM to integer conversion routines
!                    when no valid result can be expressed as a one word
!                    integer.  KFLAG = -4 is also set.

   INTEGER, PRIVATE, PARAMETER :: IUNKNO = -MXBASE*MXBASE

!             JFORM1 indicates the format used by FMOUT.

   INTEGER :: JFORM1 = 1

!             JFORM2 indicates the number of digits used in FMOUT.

   INTEGER :: JFORM2 = INITIAL_NPREC

!             KRAD = 1 indicates that trig functions use radians,
!                  = 0 means use degrees.

   INTEGER :: KRAD = 1

!             KWARN = 0 indicates that no warning message is printed
!                       and execution continues when UNKNOWN or another
!                       exception is produced.
!                   = 1 means print a warning message and continue.
!                   = 2 means print a warning message and stop.

   INTEGER :: KWARN = 1

!             KROUND = 1   Causes all results to be rounded to the
!                          nearest FM number, or to the value with
!                          an even last digit if the result is halfway
!                          between two FM numbers.
!                    = 0   Causes all results to be chopped.

   INTEGER :: KROUND = 1

CONTAINS

SUBROUTINE XP_GETS_INT (MA, I)
   TYPE (XP_REAL), INTENT (OUT) :: MA
   INTEGER, INTENT (IN) :: I
   MA = XP (I)
END SUBROUTINE

SUBROUTINE XP_GETS_REAL (MA, R)
   TYPE (XP_REAL), INTENT (OUT) :: MA
   REAL, INTENT (IN) :: R
   MA = XP (R)
END SUBROUTINE

SUBROUTINE XP_GETS_DP (MA, D)
   TYPE (XP_REAL), INTENT (OUT) :: MA
   DOUBLE PRECISION, INTENT (IN) :: D
   MA = XP (D)
END SUBROUTINE

FUNCTION XP_ABS (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_ABS
   CALL FMABS(MA%DIGITS,XP_ABS%DIGITS) ! XP_ABS = ABS(MA)
END FUNCTION

FUNCTION XP_ACOS (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_ACOS
   CALL FMACOS(MA%DIGITS,XP_ACOS%DIGITS) ! XP_COS = ACOS(MA)
END FUNCTION

FUNCTION XP_ADD (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_ADD
   CALL FMADD(MA%DIGITS,MB%DIGITS,XP_ADD%DIGITS) ! XP_ADD = MA + MB
END FUNCTION

FUNCTION XPI_ADD (MA, I)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: XPI_ADD
   XPI_ADD = MA + XP (I)
END FUNCTION

FUNCTION IXP_ADD (I, MA)
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: IXP_ADD
   IXP_ADD = MA + XP (I)
END FUNCTION

FUNCTION XP_AINT (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_AINT
   CALL FMINT(MA%DIGITS,XP_AINT%DIGITS) ! XP_AINT = AINT(MA) integer part of MA.
END FUNCTION

FUNCTION XP_ANINT (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_ANINT
   CALL FMNINT(MA%DIGITS,XP_ANINT%DIGITS) ! XP_ANINT = ANINT(MA)
END FUNCTION

FUNCTION XP_ASIN (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_ASIN
   CALL FMASIN(MA%DIGITS,XP_ASIN%DIGITS) ! XP_ASIN = ASIN(MA)
END FUNCTION

FUNCTION XP_ATAN (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_ATAN
   CALL FMATAN(MA%DIGITS,XP_ATAN%DIGITS) ! XP_ATAN = ATAN(MA)
END FUNCTION

FUNCTION XP_ATAN2 (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_ATAN2
   CALL FMATN2(MA%DIGITS,MB%DIGITS,XP_ATAN2%DIGITS) ! XP_ATAN2 = ATAN2(MA,MB)
END FUNCTION

FUNCTION XP_HUGE (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_HUGE
   CALL FMBIG(XP_HUGE%DIGITS) ! XP_HUGE = Biggest FM number less than overflow.
END FUNCTION

FUNCTION XP_EQ (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_EQ
   XP_EQ = FMCOMP (MA%DIGITS, "EQ", MB%DIGITS)
END FUNCTION

FUNCTION XP_NE (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_NE
   XP_NE = FMCOMP (MA%DIGITS, "NE", MB%DIGITS)
END FUNCTION

FUNCTION XP_LE (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_LE
   XP_LE = FMCOMP (MA%DIGITS, "LE", MB%DIGITS)
END FUNCTION

FUNCTION XP_GE (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_GE
   XP_GE = FMCOMP (MA%DIGITS, "GE", MB%DIGITS)
END FUNCTION

FUNCTION XP_LT (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_LT
   XP_LT = FMCOMP (MA%DIGITS, "LT", MB%DIGITS)
END FUNCTION

FUNCTION XP_GT (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   LOGICAL :: XP_GT
   XP_GT = FMCOMP (MA%DIGITS, "GT", MB%DIGITS)
END FUNCTION

FUNCTION XP_COS (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_COS
   CALL FMCOS(MA%DIGITS,XP_COS%DIGITS) ! XP_COS = COS(MA)
END FUNCTION

FUNCTION XP_COSH (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_COSH
   CALL FMCOSH(MA%DIGITS,XP_COSH%DIGITS) ! XP_COSH = COSH(MA)
END FUNCTION

FUNCTION XP_DIM (MA,MB)
   TYPE (XP_REAL), INTENT (IN) :: MA,MB
   TYPE (XP_REAL) :: XP_DIM
   CALL FMDIM(MA%DIGITS,MB%DIGITS,XP_DIM%DIGITS) ! XP_DIM = DIM(MA,MB)
END FUNCTION

FUNCTION XPI_DIV (MA, I)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: XPI_DIV
   CALL FMDIVI(MA%DIGITS,I,XPI_DIV%DIGITS) ! XPI_DIV = MA/I
END FUNCTION

FUNCTION IXP_DIV (I, MA)
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: IXP_DIV
   IXP_DIV = XP (I) / MA
END FUNCTION

FUNCTION XP_DIV (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_DIV
   CALL FMDIV(MA%DIGITS,MB%DIGITS,XP_DIV%DIGITS) ! XP_DIV = MA/MB
END FUNCTION

FUNCTION XP_I2M (I)
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: XP_I2M
   CALL FMI2M(I,XP_I2M%DIGITS) ! XP_I2M = X
END FUNCTION

FUNCTION XP_R2M (X)
   REAL, INTENT (IN) :: X
   TYPE (XP_REAL) :: XP_R2M
   CALL FMSP2M(X,XP_R2M%DIGITS) ! XP_R2M = X
END FUNCTION

FUNCTION XP_DP2M (D)
   DOUBLE PRECISION, INTENT (IN) :: D
   TYPE (XP_REAL) :: XP_DP2M
   CALL FMDP2M(D,XP_DP2M%DIGITS) ! XP_DP2M = X
END FUNCTION

FUNCTION XP_C2M (C)
   CHARACTER, DIMENSION (:), INTENT (IN) :: C
   TYPE (XP_REAL) :: XP_C2M
   CALL FMINP(C,XP_C2M%DIGITS,1,SIZE(C)) ! XP_C2M = LINE
        ! input conversion of C(LA) through C(LB)
        ! from characters to FM.
END FUNCTION

FUNCTION XP_EXP (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_EXP
   CALL FMEXP(MA%DIGITS,XP_EXP%DIGITS) ! XP_EXP = EXP(MA)
END FUNCTION

FUNCTION IXP_PWR (I, MA)
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: IXP_PWR
   IXP_PWR = XP (I) ** MA
END FUNCTION

FUNCTION XPI_PWR (MA,I)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: XPI_PWR
   CALL FMIPWR(MA%DIGITS,I,XPI_PWR%DIGITS) ! XPI_PWR = MA**I
END FUNCTION

FUNCTION XP_PWR (MA,MB)
   TYPE (XP_REAL), INTENT (IN) :: MA,MB
   TYPE (XP_REAL) :: XP_PWR
   CALL FMPWR(MA%DIGITS,MB%DIGITS,XP_PWR%DIGITS) ! XP_PWR = MA**MB
END FUNCTION

FUNCTION XP_LOG10(MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_LOG10
   CALL FMLG10(MA%DIGITS,XP_LOG10%DIGITS) ! XP_LOG10 = LOG10(MA)
END FUNCTION

FUNCTION XP_LOG (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_LOG
   CALL FMLN(MA%DIGITS,XP_LOG%DIGITS) ! XP_LOG = LOG(MA)
END FUNCTION

FUNCTION XP_M2DP (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   DOUBLE PRECISION  :: XP_M2DP
   CALL FMM2DP(MA%DIGITS,XP_M2DP) ! XP_M2DP  = DBLE(MA)
END FUNCTION

FUNCTION XP_M2I (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER :: XP_M2I
   CALL FMM2I(MA%DIGITS,XP_M2I) ! XP_M2I = INT(MA)
END FUNCTION

FUNCTION XP_M2SP (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   REAL :: XP_M2SP
   CALL FMM2SP(MA%DIGITS,XP_M2SP) ! XP_M2SP = REAL(MA)
END FUNCTION

FUNCTION XP_MAX (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_MAX
   CALL FMMAX(MA%DIGITS,MB%DIGITS,XP_MAX%DIGITS) ! XP_MAX = MAX(MA,MB)
END FUNCTION

FUNCTION XP_MIN (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_MIN
   CALL FMMIN(MA%DIGITS,MB%DIGITS,XP_MIN%DIGITS) ! XP_MIN = MIN(MA,MB)
END FUNCTION

FUNCTION XP_MOD (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_MOD
   CALL FMMOD(MA%DIGITS,MB%DIGITS,XP_MOD%DIGITS) ! XP_MOD = MA mod MB
END FUNCTION

FUNCTION XP_MPY (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_MPY
   CALL FMMPY(MA%DIGITS,MB%DIGITS,XP_MPY%DIGITS) ! XP_MPY = MA*MB
END FUNCTION

FUNCTION XPI_MPY (MA, I)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: XPI_MPY
   CALL FMMPYI(MA%DIGITS,I,XPI_MPY%DIGITS) ! XPI_MPY = MA*I
END FUNCTION

FUNCTION IXP_MPY (I, MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   INTEGER, INTENT (IN) :: I
   TYPE (XP_REAL) :: IXP_MPY
   CALL FMMPYI(MA%DIGITS,I,IXP_MPY%DIGITS) ! IXP_MPY = I*MA
END FUNCTION

FUNCTION XP_PI (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_PI
   CALL FMPI(XP_PI%DIGITS) ! XP_PI = pi
END FUNCTION

FUNCTION XP_SIGN (MA,MB)
   TYPE (XP_REAL), INTENT (IN) :: MA,MB
   TYPE (XP_REAL) :: XP_SIGN
   CALL FMSIGN(MA%DIGITS,MB%DIGITS,XP_SIGN%DIGITS) ! XP_SIGN = SIGN(MA,MB)
END FUNCTION

FUNCTION XP_SIN (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_SIN
   CALL FMSIN(MA%DIGITS,XP_SIN%DIGITS) ! XP_SIN = SIN(MA)
END FUNCTION

FUNCTION XP_SINH (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_SINH
   CALL FMSINH(MA%DIGITS,XP_SINH%DIGITS) ! XP_SING = SINH(MA)
END FUNCTION

FUNCTION XP_SQRT (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_SQRT
   CALL FMSQRT(MA%DIGITS,XP_SQRT%DIGITS) ! XP_SQRT = SQRT(MA)
END FUNCTION

FUNCTION XP_SUB (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_SUB
   CALL FMSUB(MA%DIGITS,MB%DIGITS,XP_SUB%DIGITS) ! XP_SUB = MA - MB
END FUNCTION

FUNCTION XP_TAN (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_TAN
   CALL FMTAN(MA%DIGITS,XP_TAN%DIGITS) ! XP_TAN = TAN(MA)
END FUNCTION

FUNCTION XP_TANH (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   TYPE (XP_REAL) :: XP_TANH
   CALL FMTANH(MA%DIGITS,XP_TANH%DIGITS) ! XP_TANH = TANH(MA)
END FUNCTION

SUBROUTINE XP_PRINT (MA)
   TYPE (XP_REAL), INTENT (IN) :: MA
   CALL FMPRNT (MA%DIGITS)
END SUBROUTINE

SUBROUTINE XP_SET (NPREC)

! Sets number of digits of precision
      NPSAVE = NPREC
      NDIG = 2 + (NPREC+2)/K_RANGE
      IF (NDIG.LT.2 .OR. NDIG.GT.MXNDIG) THEN
          NDIG = MAX(2,MIN(MXNDIG,NDIG))
          WRITE (KW,120) NPREC,NDIG
  120     FORMAT(//' PRECISION OUT OF RANGE WHEN CALLING XP_SET.',  &
                 '  NPREC =',I20/' THE NEAREST VALID NDIG WILL BE USED',  &
                 ' INSTEAD:   NDIG =',I6//)
          NPSAVE = 0
      ENDIF

END SUBROUTINE

SUBROUTINE FMABS(MA,MB)

!  MB = ABS(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 1
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)

      KFLAG = 0
      KWSV = KWARN
      KWARN = 0
      CALL FMEQU(MA,MB,NDIG,NDIG)
      MB(2) = ABS(MB(2))
      KWARN = KWSV

      IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMACOS(MA,MB)

!  MB = ACOS(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(2,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      MA2 = MA(2)
      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Use ACOS(X) = ATAN(SQRT(1-X*X)/X)

      MB(2) = ABS(MB(2))
      CALL FMI2M(1,M05)
      CALL FMSUB(M05,MB,M03)
      CALL FMADD(M05,MB,M04)
      CALL FMMPY(M03,M04,M04)
      CALL FMSQRT(M04,M04)
      CALL FMDIV(M04,MB,MB)

      CALL FMATAN(MB,MB)

      IF (MA2.LT.0) THEN
          IF (KRAD.EQ.1) THEN
              CALL FMPI(M05)
          ELSE
              CALL FMI2M(180,M05)
          ENDIF
          CALL FMSUB(M05,MB,MB)
      ENDIF

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMADD(MA,MB,MC)

!  MC = MA + MB

!  This routine performs the trace printing for addition.
!  FMADD2 is used to do the arithmetic.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 3
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      CALL FMADD2(MA,MB,MC)

      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMADD2(MA,MB,MC)

!  Internal addition routine.  MC = MA + MB

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      IF (KARGSW.NE.1) THEN
          CALL FMARGS(3,2,MA,MB,KRESLT)
          IF (KRESLT.NE.0) THEN
              CALL FMRSLT(MA,MB,MC,KRESLT)
              RETURN
          ENDIF
      ELSE
          IF (MA(2).EQ.0) THEN
              CALL FMEQU(MB,MC,NDIG,NDIG)
              KFLAG = 1
              RETURN
          ENDIF
          IF (MB(2).EQ.0) THEN
              CALL FMEQU(MA,MC,NDIG,NDIG)
              KFLAG = 1
              RETURN
          ENDIF
      ENDIF

      KFLAG = 0
      N1 = NDIG + 1

!             NGUARD is the number of guard digits used.

      IF (NCALL.GT.1) THEN
          NGUARD = 1
      ELSE
          B = JBASE
          NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0
          IF (NGUARD.GT.NDIG) NGUARD = NDIG
      ENDIF
      NMWA = N1 + NGUARD

!             Save the signs of MA and MB and then work with
!             positive numbers.
!             JSIGN is the sign of the result of MA + MB.

      JSIGN = 1
      MA2 = MA(2)
      MB2 = MB(2)
      MA(2) = ABS(MA(2))
      MB(2) = ABS(MB(2))

!             See which one is larger in absolute value.

      IF (MA(1).GT.MB(1)) THEN
          JCOMP = 1
          GO TO 120
      ENDIF
      IF (MB(1).GT.MA(1)) THEN
          JCOMP = 3
          GO TO 120
      ENDIF
      NLAST = NDIG + 1

      DO 110 J = 2, NLAST
         IF (ABS(MA(J)).GT.ABS(MB(J))) THEN
             JCOMP = 1
             GO TO 120
         ENDIF
         IF (ABS(MB(J)).GT.ABS(MA(J))) THEN
             JCOMP = 3
             GO TO 120
         ENDIF
  110 CONTINUE

      JCOMP = 2

  120 IF (JCOMP.LT.3) THEN
          KBIGMA = 1
          IF (MA2.LT.0) JSIGN = -1
          IF (MA2*MB2.GT.0) THEN
              CALL FMADDP(MA,MB,NMWA)
          ELSE
              CALL FMADDN(MA,MB,NGUARD,NMWA)
          ENDIF
      ELSE
          KBIGMA = 0
          IF (MB2.LT.0) JSIGN = -1
          IF (MA2*MB2.GT.0) THEN
              CALL FMADDP(MB,MA,NMWA)
          ELSE
              CALL FMADDN(MB,MA,NGUARD,NMWA)
          ENDIF
      ENDIF
      MA(2) = MA2
      MB(2) = MB2


!             Round the result.

      CALL FMRND(NDIG,NGUARD,0)

!             See if the result is equal to one of the input arguments.

      K = ABS(MA(1)-MB(1))
      IF (K.LT.NDIG) GO TO 150
      IF (K.GT.NDIG+1) THEN
          KFLAG = 1
          GO TO 150
      ENDIF

      N2 = NDIG + 4
      IF (KBIGMA.EQ.1) THEN
          DO 130 J = 3, N1
             IF (MWA(N2-J).NE.MA(N2-J)) GO TO 150
  130     CONTINUE
          IF (MWA(1).NE.MA(1)) GO TO 150
          IF (MWA(2).NE.ABS(MA(2))) GO TO 150
      ELSE
          DO 140 J = 3, N1
             IF (MWA(N2-J).NE.MB(N2-J)) GO TO 150
  140     CONTINUE
          IF (MWA(1).NE.MB(1)) GO TO 150
          IF (MWA(2).NE.ABS(MB(2))) GO TO 150
      ENDIF
      KFLAG = 1

!             Transfer to MC and fix the sign of the result.

  150 CALL FMMOVE(MC)
      IF (JSIGN.LT.0) MC(2) = -MC(2)

END SUBROUTINE

SUBROUTINE FMADDN(MA,MB,NGUARD,NMWA)

!  Internal addition routine.  MWA = MA - MB
!  The arguments are such that MA.GE.MB.GE.0.

!  NGUARD is the number of guard digits being carried.
!  NMWA is the number of words in MWA which will be used.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      N1 = NDIG + 1

!             Check for an insignificant operand.

      K = MA(1) - MB(1)
      IF (K.GE.NDIG+2) THEN
          DO 110 J = 1, N1
             MWA(J) = MA(J)
  110     CONTINUE
          MWA(N1+1) = 0
          RETURN
      ENDIF
      IF (NGUARD.LE.1) NMWA = N1 + 2

!             Subtract MB from MA.

      KP1 = MIN(N1,K+1)
      MWA(K+1) = 0
      DO 120 J = 1, KP1
         MWA(J) = MA(J)
  120 CONTINUE
      KP2 = K + 2
      IF (KP2.LE.N1) THEN
          DO 130 J = KP2, N1
             MWA(J) = MA(J) - MB(J-K)
  130     CONTINUE
      ENDIF
      N2 = NDIG + 2
      IF (N2-K.LE.1) N2 = 2 + K
      NK = MIN(NMWA,N1+K)
      IF (N2.LE.NK) THEN
          DO 140 J = N2, NK
             MWA(J) = -MB(J-K)
  140     CONTINUE
      ENDIF
      NK1 = NK + 1
      IF (NK1.LE.NMWA) THEN
          DO 150 J = NK1, NMWA
             MWA(J) = 0
  150     CONTINUE
      ENDIF

!             Normalize.  Fix the sign of any negative digit.

      IF (K.GT.0) THEN
          KB = NMWA - KP2 + 1
          K2 = NMWA + 1
          DO 160 J = 1, KB
             IF (MWA(K2-J).LT.0) THEN
                 MWA(K2-J) = MWA(K2-J) + JBASE
                 MWA(NMWA-J) = MWA(NMWA-J) - 1
             ENDIF
  160     CONTINUE
          KPT = KP2
  170     KPT = KPT - 1
          IF (MWA(KPT).LT.0 .AND. KPT.GE.3) THEN
              MWA(KPT) = MWA(KPT) + JBASE
              MWA(KPT-1) = MWA(KPT-1) - 1
              GO TO 170
          ENDIF
          GO TO 190
      ENDIF

      IF (K.EQ.0) THEN
          KP = N1 + 3
          KPM1 = KP - 1
          DO 180 J = 3, N1
             IF (MWA(KP-J).LT.0) THEN
                 MWA(KP-J) = MWA(KP-J) + JBASE
                 MWA(KPM1-J) = MWA(KPM1-J) - 1
             ENDIF
  180     CONTINUE
      ENDIF

!             Shift left if there are any leading zeros in the mantissa.

  190 DO 200 J = 2, NMWA
         KSH = J - 2
         IF (MWA(J).GT.0) GO TO 210
  200 CONTINUE
      MWA(1) = 0
      RETURN

  210 KL = NMWA - KSH
      IF (KSH.GT.0) THEN
          DO 220 J = 2, KL
             L = J + KSH
             MWA(J) = MWA(L)
  220     CONTINUE
          KL = KL + 1
          DO 230 J = KL, NMWA
             MWA(J) = 0
  230     CONTINUE
          MWA(1) = MWA(1) - KSH
      ENDIF

END SUBROUTINE

SUBROUTINE FMADDP(MA,MB,NMWA)

!  Internal addition routine.  MWA = MA + MB
!  The arguments are such that MA.GE.MB.GE.0.

!  NMWA is the number of words in MWA which will be used.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      N1 = NDIG + 1

!             Check for an insignificant operand.

      K = MA(1) - MB(1)
      IF (K.GE.NDIG+1) THEN
          DO 110 J = 1, N1
             MWA(J) = MA(J)
  110     CONTINUE
          MWA(N1+1) = 0
          RETURN
      ENDIF

!             Add MA and MB.

      KP1 = K + 1
      DO 120 J = 1, KP1
         MWA(J) = MA(J)
  120 CONTINUE
      KP2 = K + 2
      IF (KP2.LE.N1) THEN
          DO 130 J = KP2, N1
             MWA(J) = MA(J) + MB(J-K)
  130     CONTINUE
      ENDIF
      N2 = NDIG + 2
      NK = MIN(NMWA,N1+K)
      IF (N2.LE.NK) THEN
          DO 140 J = N2, NK
             MWA(J) = MB(J-K)
  140     CONTINUE
      ENDIF
      NK1 = NK + 1
      IF (NK1.LE.NMWA) THEN
          DO 150 J = NK1, NMWA
             MWA(J) = 0
  150     CONTINUE
      ENDIF

!             Normalize.  Fix any digit not less than JBASE.

      IF (K.EQ.NDIG) RETURN

      IF (K.GT.0) THEN
          KB = N1 - KP2 + 1
          N2 = N1 + 1
          DO 160 J = 1, KB
             IF (MWA(N2-J).GE.JBASE) THEN
                 MWA(N2-J) = MWA(N2-J) - JBASE
                 MWA(N1-J) = MWA(N1-J) + 1
             ENDIF
  160     CONTINUE
          KPT = KP2
  170     KPT = KPT - 1
          IF (MWA(KPT).GE.JBASE .AND. KPT.GE.3) THEN
              MWA(KPT) = MWA(KPT) - JBASE
              MWA(KPT-1) = MWA(KPT-1) + 1
              GO TO 170
          ENDIF
          GO TO 190
      ENDIF

      IF (K.EQ.0) THEN
          KP = N1 + 3
          KPM1 = KP - 1
          DO 180 J = 3, N1
             IF (MWA(KP-J).GE.JBASE) THEN
                 MWA(KP-J) = MWA(KP-J) - JBASE
                 MWA(KPM1-J) = MWA(KPM1-J) + 1
             ENDIF
  180     CONTINUE
      ENDIF

!             Shift right if the leading digit is not less than JBASE.

  190 IF (MWA(2).GE.JBASE) THEN
  200     KP = NMWA + 4
          DO 210 J = 4, NMWA
             MWA(KP-J) = MWA(KP-J-1)
  210     CONTINUE
          KT = MWA(2)/JBASE
          MWA(3) = MWA(2) - KT*JBASE
          MWA(2) = KT
          MWA(1) = MWA(1) + 1
          IF (MWA(2).GE.JBASE) GO TO 200
      ENDIF

END SUBROUTINE

SUBROUTINE FMARGS(KROUTN,NARGS,MA,MB,KRESLT)

!  Check the input arguments to a routine for special cases.

!  KROUTN - ID number of the subroutine which was called
!  NARGS  - The number of input arguments (1 or 2)
!  MA     - First input argument
!  MB     - Second input argument (if NARGS is 2)
!  KRESLT - Result code returned to the calling routine.

!  Result codes:

!   0 - Perform the normal operation
!   1 - The result is the first input argument
!   2 - The result is the second input argument
!   3 - The result is -OVERFLOW
!   4 - The result is +OVERFLOW
!   5 - The result is -UNDERFLOW
!   6 - The result is +UNDERFLOW
!   7 - The result is -1.0
!   8 - The result is +1.0
!   9 - The result is -pi/2
!  10 - The result is +pi/2
!  11 - The result is 0.0
!  12 - The result is UNKNOWN
!  13 - The result is +pi
!  14 - The result is -pi/4
!  15 - The result is +pi/4

      DIMENSION MA(LUNPCK),MB(LUNPCK),KADD(15,15),KMPY(15,15),  &
                KDIV(15,15),KPWR(15,15),KSQRT(15),KEXP(15),KLN(15),  &
                KSIN(15),KCOS(15),KTAN(15),KASIN(15),KACOS(15),  &
                KATAN(15),KSINH(15),KCOSH(15),KTANH(15),KLG10(15)

!             These tables define the result codes to be returned for
!             given values of the input argument(s).

!             For example, in row 7 column 2 of this DATA statement
!             KADD(2,7) = 2 means that if the first argument in a call
!             to FMADD is in category 7 ( -UNDERFLOW ) and the second
!             argument is in category 2 ( near -OVERFLOW but
!             representable ) then the result code is 2 ( the value
!             of the sum is equal to the second input argument).
!             See routine FMCAT for descriptions of the categories.

      DATA KADD/ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,12,12,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0,12,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0,12, 1,12, 0, 0, 0, 0, 0, 4,  &
                 3, 2, 2, 2, 2,12,12, 5,12,12, 2, 2, 2, 2, 4,  &
                 3, 2, 2, 2, 2, 2, 5, 2, 6, 2, 2, 2, 2, 2, 4,  &
                 3, 2, 2, 2, 2,12,12, 6,12,12, 2, 2, 2, 2, 4,  &
                 3, 0, 0, 0, 0, 0,12, 1,12, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                12, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 4,  &
                12,12, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 /

      DATA KMPY/ 4, 4, 4, 4,12,12,12,11,12,12,12, 3, 3, 3, 3,  &
                 4, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 3,  &
                 4, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 3,  &
                 4, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0, 3,  &
                12, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0,12,  &
                12, 0, 0, 0, 0, 0, 6,11, 5, 0, 0, 1, 0, 0,12,  &
                12,12,12, 6, 6, 6, 6,11, 5, 5, 5, 5,12,12,12,  &
                11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,  &
                12,12,12, 5, 5, 5, 5,11, 6, 6, 6, 6,12,12,12,  &
                12, 0, 0, 0, 0, 0, 5,11, 6, 0, 0, 1, 0, 0,12,  &
                12, 0, 0, 0, 0, 0, 5,11, 6, 0, 0, 1, 0, 0,12,  &
                 3, 2, 2, 2, 2, 2, 5,11, 6, 2, 2, 2, 2, 2, 4,  &
                 3, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 4,  &
                 3, 0, 0, 0, 0, 0,12,11,12, 0, 0, 1, 0, 0, 4,  &
                 3, 3, 3, 3,12,12,12,11,12,12,12, 4, 4, 4, 4 /

      DATA KDIV/12,12,12, 4, 4, 4, 4,12, 3, 3, 3, 3,12,12,12,  &
                12, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0,12,  &
                12, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0,12,  &
                 6, 0, 0, 0, 0, 0, 4,12, 3, 0, 0, 1, 0, 0, 5,  &
                 6, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 5,  &
                 6, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 5,  &
                 6, 6, 6, 6,12,12,12,12,12,12,12, 5, 5, 5, 5,  &
                11,11,11,11,11,11,11,12,11,11,11,11,11,11,11,  &
                 5, 5, 5, 5,12,12,12,12,12,12,12, 6, 6, 6, 6,  &
                 5, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 6,  &
                 5, 0, 0, 0, 0, 0,12,12,12, 0, 0, 1, 0, 0, 6,  &
                 5, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0, 6,  &
                12, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0,12,  &
                12, 0, 0, 0, 0, 0, 3,12, 4, 0, 0, 1, 0, 0,12,  &
                12,12,12, 3, 3, 3, 3,12, 4, 4, 4, 4,12,12,12 /

      DATA KPWR/12,12, 0, 5,12,12,12, 8,12,12,12, 3, 0,12,12,  &
                12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12,  &
                12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12,  &
                12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12,  &
                12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12,  &
                12,12, 0, 0,12,12,12, 8,12,12,12, 1, 0,12,12,  &
                12,12, 0, 3,12,12,12, 8,12,12,12, 5, 0,12,12,  &
                12,12,12,12,12,12,12,12,11,11,11,11,11,11,11,  &
                 4, 4, 4, 4,12,12,12, 8,12,12,12, 6, 6, 6, 6,  &
                 4, 4, 0, 0, 0, 8, 8, 8, 8, 0, 0, 1, 0, 6, 6,  &
                 4, 4, 0, 0, 0, 8, 8, 8, 8, 0, 0, 1, 0, 6, 6,  &
                 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,  &
                 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 1, 0, 4, 4,  &
                 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 1, 0, 4, 4,  &
                 6, 6, 6, 6,12,12,12, 8,12,12,12, 4, 4, 4, 4 /

      DATA KSQRT/12,12,12,12,12,12,12,11,12, 0, 0, 8, 0, 0,12/
      DATA KEXP / 6, 6, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0, 4, 4/
      DATA KLN  /12,12,12,12,12,12,12,12,12, 0, 0,11, 0, 0,12/
      DATA KSIN /12,12, 0, 0, 0, 0, 5,11, 6, 0, 0, 0, 0,12,12/
      DATA KCOS /12,12, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0,12,12/
      DATA KTAN /12,12, 0, 0, 0, 0, 5,11, 6, 0, 0, 0, 0,12,12/
      DATA KASIN/12,12,12, 9, 0, 0, 5,11, 6, 0, 0,10,12,12,12/
      DATA KACOS/12,12,12,13, 0,10,10,10,10,10, 0,11,12,12,12/
      DATA KATAN/ 9, 9, 0,14, 0, 0, 5,11, 6, 0, 0,15, 0,10,10/
      DATA KSINH/ 3, 3, 0, 0, 0, 1, 5,11, 6, 1, 0, 0, 0, 4, 4/
      DATA KCOSH/ 4, 4, 0, 0, 0, 8, 8, 8, 8, 8, 0, 0, 0, 4, 4/
      DATA KTANH/ 7, 7, 0, 0, 0, 1, 5,11, 6, 1, 0, 0, 0, 8, 8/
      DATA KLG10/12,12,12,12,12,12,12,12,12, 0, 0,11, 0, 0,12/

      KRESLT = 12
      KFLAG = -4
      IF (MA(1).EQ.KUNKNO) RETURN
      IF (NARGS.EQ.2) THEN
          IF (MB(1).EQ.KUNKNO) RETURN
      ENDIF
      KFLAG = 0

!             Check the validity of parameters if this is a user call.

      IF (NCALL.GT.1) GO TO 130

!             Check NDIG.

      IF (NDIG.LT.2 .OR. NDIG.GT.MXNDIG) THEN
          KFLAG = -1
          CALL FMWARN
          NDS = NDIG
          IF (NDIG.LT.2) NDIG = 2
          IF (NDIG.GT.MXNDIG) NDIG = MXNDIG
          WRITE (KW,110) NDS,NDIG
  110     FORMAT(' NDIG WAS',I10,'.  IT HAS BEEN CHANGED TO',I10,'.')
          RETURN
      ENDIF

!             Check JBASE.

      IF (JBASE.LT.2 .OR. JBASE.GT.MXBASE) THEN
          KFLAG = -2
          CALL FMWARN
          JBS = JBASE
          IF (JBASE.LT.2) JBASE = 2
          IF (JBASE.GT.MXBASE) JBASE = MXBASE
          WRITE (KW,120) JBS,JBASE
  120     FORMAT(' JBASE WAS',I10,'.  IT HAS BEEN CHANGED TO',I10,'.')
          RETURN
      ENDIF

!             Check exponent range.

      IF (MA(1).GT.MXEXP+1 .OR. MA(1).LT.-MXEXP) THEN
          IF (ABS(MA(1)).NE.KEXPOV .OR. ABS(MA(2)).NE.1) THEN
              CALL FMIM(0,MA)
              KFLAG = -3
              CALL FMWARN
              MA(1) = KUNKNO
              MA(2) = 1
              RETURN
          ENDIF
      ENDIF
      IF (NARGS.EQ.2) THEN
          IF (MB(1).GT.MXEXP+1 .OR. MB(1).LT.-MXEXP) THEN
              IF (ABS(MB(1)).NE.KEXPOV .OR. ABS(MB(2)).NE.1) THEN
                  CALL FMIM(0,MB)
                  KFLAG = -3
                  CALL FMWARN
                  MB(1) = KUNKNO
                  MB(2) = 1
                  RETURN
              ENDIF
          ENDIF
      ENDIF

!             Check for special cases.

  130 CALL FMCAT(MA,NCATMA)
      NCATMB = 0
      IF (NARGS.EQ.2) CALL FMCAT(MB,NCATMB)

      SELECT CASE (KROUTN)
      CASE (3)
          KRESLT = KADD(NCATMB,NCATMA)

      CASE (29)
          KRESLT = KMPY(NCATMB,NCATMA)

      CASE (12)
          KRESLT = KDIV(NCATMB,NCATMA)

      CASE (34)
          KRESLT = KPWR(NCATMB,NCATMA)

      CASE (39)
          KRESLT = KSQRT(NCATMA)

      CASE (15)
          KRESLT = KEXP(NCATMA)

      CASE (21)
          KRESLT = KLN(NCATMA)

      CASE (36)
          KRESLT = KSIN(NCATMA)

      CASE (9)
          KRESLT = KCOS(NCATMA)

      CASE (41)
          KRESLT = KTAN(NCATMA)

      CASE (4)
          KRESLT = KASIN(NCATMA)
          IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12

      CASE (2)
          KRESLT = KACOS(NCATMA)

      CASE (5)
          KRESLT = KATAN(NCATMA)
          IF ((NCATMA.EQ.7.OR.NCATMA.EQ.9) .AND. KRAD.EQ.0) KRESLT = 12

      CASE (37)
          KRESLT = KSINH(NCATMA)

      CASE (10)
          KRESLT = KCOSH(NCATMA)

      CASE (42)
          KRESLT = KTANH(NCATMA)

      CASE (20)
          GO TO 140

      CASE DEFAULT
      KRESLT = 0

      END SELECT

      RETURN

  140 IF (KRESLT.EQ.12) THEN
          KFLAG = -4
          CALL FMWARN
      ENDIF
      IF (KRESLT.EQ.3 .OR. KRESLT.EQ.4) THEN
          IF (NCATMA.EQ.1 .OR. NCATMA.EQ.7 .OR. NCATMA.EQ.9 .OR.  &
              NCATMA.EQ.15 .OR. NCATMB.EQ.1 .OR. NCATMB.EQ.7 .OR.  &
              NCATMB.EQ.9 .OR. NCATMB.EQ.15) THEN
              KFLAG = -5
          ELSE
              KFLAG = -5
              CALL FMWARN
          ENDIF
      ENDIF
      IF (KRESLT.EQ.5 .OR. KRESLT.EQ.6) THEN
          IF (NCATMA.EQ.1 .OR. NCATMA.EQ.7 .OR. NCATMA.EQ.9 .OR.  &
              NCATMA.EQ.15 .OR. NCATMB.EQ.1 .OR. NCATMB.EQ.7 .OR.  &
              NCATMB.EQ.9 .OR. NCATMB.EQ.15) THEN
              KFLAG = -6
          ELSE
              KFLAG = -6
              CALL FMWARN
          ENDIF
      ENDIF
END SUBROUTINE

SUBROUTINE FMASIN(MA,MB)

!  MB = ARCSIN(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(4,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Use ASIN(X) = ATAN(X/SQRT(1-X*X))

      CALL FMI2M(1,M05)
      CALL FMSUB(M05,MB,M03)
      CALL FMADD(M05,MB,M04)
      CALL FMMPY(M03,M04,M04)
      CALL FMSQRT(M04,M04)
      CALL FMDIV(MB,M04,MB)

      CALL FMATAN(MB,MB)

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMATAN(MA,MB)

!  MB = ARCTAN(MA)

      DOUBLE PRECISION X,XB,XM
      DIMENSION MA(LUNPCK),MB(LUNPCK),NSTACK(19)

      CALL FMENTR(5,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,M05,NDSAVE,NDIG)

!             If MA.GE.1 work with 1/MA.

      MA1 = MA(1)
      MA2 = MA(2)
      M05(2) = ABS(M05(2))
      IF (MA1.GE.1) THEN
          CALL FMI2M(1,MB)
          CALL FMDIV(MB,M05,M05)
      ENDIF

      KRSAVE = KRAD
      KRAD = 1
      KWSV = KWARN

      X = M05(1)
      XB = JBASE
      XM = MXBASE

!             In case pi has not been computed at the current precision
!             and will be needed here, get it to full precision first
!             to avoid repeated calls at increasing precision during
!             Newton iteration.

      IF (MA1.GE.1 .OR. KRSAVE.EQ.0) THEN
          IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMPI2(MPISAV)
              NJBPI = JBASE
              NDIGPI = NDIG
              NDIG = NDSV
          ENDIF
      ENDIF

!             If the argument is small, use the Taylor series,
!             otherwise use Newton iteration.

      IF (X*LOG(XB).LT.-5.0*LOG(XM)) THEN
          KWARN = 0
          CALL FMEQU(M05,MB,NDIG,NDIG)
          CALL FMMPY(M05,M05,M06)
          J = 3
          NDSAV1 = NDIG

  110     CALL FMMPY(M05,M06,M05)
          M05(2) = -M05(2)
          CALL FMDIVI(M05,J,M03)
          NDIG = NDSAV1
          CALL FMADD(MB,M03,MB)
          IF (KFLAG.EQ.1) THEN
              KFLAG = 0
              GO TO 130
          ENDIF
          NDIG = NDSAV1 - (MB(1)-M03(1))
          IF (NDIG.LT.2) NDIG = 2
          J = J + 2
          GO TO 110
      ELSE

          CALL FMM2DP(M05,X)
          X = ATAN(X)
          CALL FMDP2M(X,MB)
          CALL FMDIG(NSTACK,KST)

!             Newton iteration.

          DO 120 J = 1, KST
             NDIG = NSTACK(J)
             CALL FMSIN(MB,M06)
             CALL FMMPY(M06,M06,M03)
             CALL FMI2M(1,M04)
             CALL FMSUB(M04,M03,M03)
             CALL FMSQRT(M03,M04)
             CALL FMDIV(M06,M04,M04)
             CALL FMSUB(M04,M05,M04)
             CALL FMMPY(M03,M04,M04)
             CALL FMSUB(MB,M04,MB)
  120     CONTINUE
      ENDIF

!             If MA.GE.1 use pi/2 - ATAN(1/MA)

  130 IF (MA1.GE.1) THEN
          CALL FMDIVI(MPISAV,2,M06)
          CALL FMSUB(M06,MB,MB)
      ENDIF

!             Convert to degrees if necessary, round and return.

      KRAD = KRSAVE
      IF (KRAD.EQ.0) THEN
          CALL FMMPYI(MB,180,MB)
          CALL FMDIV(MB,MPISAV,MB)
      ENDIF
      IF (MA2.LT.0) MB(2) = -MB(2)

      IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMATN2(MA,MB,MC)

!  MC = ATAN2(MA,MB)

!  MC is returned as the angle between -pi and pi (or -180 and 180 if
!  degree mode is selected) for which TAN(MC) = MA/MB.  MC is an angle
!  for the point (MB,MA) in polar coordinates.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      CALL FMENTR(6,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      KARGSW = 0
      KWSV = KWARN
      KWARN = 0
      CALL FMEQU(MA,M01,NDSAVE,NDIG)
      CALL FMEQU(MB,M02,NDSAVE,NDIG)

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO .OR.  &
         (MA(2).EQ.0 .AND. MB(2).EQ.0)) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
          GO TO 110
      ENDIF

      IF (MB(2).EQ.0 .AND. MA(2).GT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(2).EQ.0 .AND. MA(2).LT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(-90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,-2,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MA(1).EQ.KEXPOV .AND. MB(1).LT.MXSAVE-NDIG-2) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

      IF (MA(1).EQ.KEXPUN .AND. (-MB(1)).LT.MXSAVE-NDIG-2 .AND.  &
                                 MB(2).LT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(180,MC)
          ELSE
              CALL FMPI(MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

      IF (MB(1).EQ.KEXPOV .AND. MA(1).LT.MXSAVE-NDIG-2 .AND.  &
                                MB(2).LT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(180,MC)
          ELSE
              CALL FMPI(MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

      IF (MB(1).EQ.KEXPUN .AND. MA(2).EQ.0) THEN
          IF (MB(2).LT.0) THEN
              IF (KRAD.EQ.0) THEN
                  CALL FMI2M(180,MC)
              ELSE
                  CALL FMPI(MC)
              ENDIF
          ELSE
              CALL FMI2M(0,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(1).EQ.KEXPUN .AND. (-MA(1)).LT.MXSAVE-NDIG-2) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

!             Determine the quadrant for the result, then use FMATAN.

      IF (MA(2).GE.0 .AND. MB(2).GT.0) JQUAD = 1
      IF (MA(2).GE.0 .AND. MB(2).LT.0) JQUAD = 2
      IF (MA(2).LT.0 .AND. MB(2).LT.0) JQUAD = 3
      IF (MA(2).LT.0 .AND. MB(2).GT.0) JQUAD = 4

      CALL FMDIV(M01,M02,MC)
      MC(2) = ABS(MC(2))
      CALL FMATAN(MC,MC)

      IF (JQUAD.EQ.2 .OR. JQUAD.EQ.3) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(180,M05)
              CALL FMSUB(M05,MC,MC)
          ELSE
              CALL FMPI(M05)
              CALL FMSUB(M05,MC,MC)
          ENDIF
      ENDIF

      IF ((JQUAD.EQ.3 .OR. JQUAD.EQ.4) .AND. MC(1).NE.KUNKNO)  &
          MC(2) = -MC(2)

!             Round the result and return.

  110 IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMBIG(MA)

!     MA = The biggest representable FM number using the current base
!          and precision.
!          The smallest positive number is then 1.0/MA.
!          Because of rounding, 1.0/(1.0/MA) will then overflow.

      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 7

      KFLAG = 0
      N1 = NDIG + 1
      DO 110 J = 2, N1
         MA(J) = JBASE - 1
  110 CONTINUE
      MA(1) = MXEXP + 1

      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMCAT(MA,NCAT)

!  NCAT is returned as the category of MA.  This is used by the various
!  arithmetic routines to handle special cases such as:
!  'number greater than 1' + 'underflowed result' is the first argument,
!  'overflowed result' / 'overflowed result' is 'unknown'.

!  NCAT       range

!   1.         -OV                OV stands for overflowed results.
!   2.   (-OV   , -OVTH)             ( MA(1) .GE. MAXEXP+2 )
!   3.   (-OVTH ,    -1)
!   4.         -1                 OVTH stands for a representable
!   5.   (-1    , -UNTH)               number near the overflow
!   6.   (-UNTH ,   -UN)               threshold.
!   7.         -UN                     ( MA(1) .GE. MAXEXP-NDIG+1 )
!   8.          0
!   9.         +UN                UN stands for underflowed results.
!  10.   (+UN   , +UNTH)             ( MA(1) .LE. -MAXEXP-1 )
!  11.   (+UNTH ,    +1)
!  12.         +1                 UNTH stands for a representable
!  13.   (+1    , +OVTH)               number near the underflow
!  14.   (+OVTH ,   +OV)               threshold.
!  15.         +OV                     ( MA(1) .LE. -MAXEXP+NDIG-1 )
!  16.       UNKNOWN

      DIMENSION MA(LUNPCK)

!             Check for special symbols.

      NCAT = 16
      IF (MA(1).EQ.KUNKNO) RETURN

      IF (MA(1).EQ.KEXPOV) THEN
          NCAT = 15
          IF (MA(2).LT.0) NCAT = 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPUN) THEN
          NCAT = 9
          IF (MA(2).LT.0) NCAT = 7
          RETURN
      ENDIF

      IF (MA(2).EQ.0) THEN
          NCAT = 8
          RETURN
      ENDIF

!             Check for +1 or -1.

      MA2 = ABS(MA(2))
      IF (MA(1).EQ.1 .AND. MA2.EQ.1) THEN
          NLAST = NDIG + 1
          IF (NLAST.GE.3) THEN
              DO 110 J = 3, NLAST
                 IF (MA(J).NE.0) GO TO 120
  110         CONTINUE
          ENDIF
          NCAT = 12
          IF (MA(2).LT.0) NCAT = 4
          RETURN
      ENDIF

  120 IF (MA(1).GE.MXEXP-NDIG+1) THEN
          NCAT = 14
          IF (MA(2).LT.0) NCAT = 2
          RETURN
      ENDIF

      IF (MA(1).GE.1) THEN
          NCAT = 13
          IF (MA(2).LT.0) NCAT = 3
          RETURN
      ENDIF

      IF (MA(1).GE.-MXEXP+NDIG) THEN
          NCAT = 11
          IF (MA(2).LT.0) NCAT = 5
          RETURN
      ENDIF

      IF (MA(1).GE.-MXEXP) THEN
          NCAT = 10
          IF (MA(2).LT.0) NCAT = 6
          RETURN
      ENDIF

END SUBROUTINE

FUNCTION FMCOMP(MA,LREL,MB)
   LOGICAL FMCOMP
   DIMENSION MA(LUNPCK),MB(LUNPCK)
   CHARACTER *2 LREL

!  Logical comparison of FM numbers MA and MB.

!  LREL is a CHARACTER *2 description of the comparison to be done:
!  LREL = 'EQ' returns FMCOMP = .TRUE. if MA.EQ.MB
!       = 'NE', 'GE', 'GT', 'LE', 'LT' also work like a logical IF.

!  For comparisons involving 'UNKNOWN' or two identical special symbols
!  such as +OVERFLOW,'EQ',+OVERFLOW FMCOMP is returned FALSE and a
!  KFLAG = -4 error condition is returned.

      NCALL = NCALL + 1
      KSTACK(NCALL) = 8

      IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.2) THEN
          WRITE (KW,110)
  110     FORMAT(' INPUT TO FMCOMP')

          IF (NTRACE.GT.0) THEN
              CALL FMPRNT(MA)
              WRITE (KW,120) LREL
  120         FORMAT(7X,'.',A2,'.')
              CALL FMPRNT(MB)
          ELSE
              CALL FMNTRJ(MA,NDIG)
              WRITE (KW,120) LREL
              CALL FMNTRJ(MB,NDIG)
          ENDIF
      ENDIF

!             JCOMP will be 1 if MA.GT.MB
!                           2 if MA.EQ.MB
!                           3 if MA.LT.MB

!             Check for special cases.

      IF (LREL.NE.'EQ' .AND. LREL.NE.'NE' .AND. LREL.NE.'LT' .AND.  &
          LREL.NE.'GT' .AND. LREL.NE.'LE' .AND. LREL.NE.'GE') THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          IF (KWARN.LE.0) GO TO 170
          WRITE (KW,130) LREL
  130     FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE',  &
                 ' FMCOMP'//1X,A,' IS NOT ONE OF THE SIX RECOGNIZED',  &
                 ' COMPARISONS.'//' .FALSE. HAS BEEN RETURNED.'/)
          GO TO 170
      ENDIF

      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          GO TO 170
      ENDIF

      IF (ABS(MA(1)).EQ.KEXPOV .AND. MA(1).EQ.MB(1) .AND.  &
          MA(2).EQ.MB(2)) THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          IF (KWARN.LE.0) GO TO 170
          WRITE (KW,140)
  140     FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE',  &
                 ' FMCOMP'//' TWO NUMBERS IN THE SAME OVERFLOW OR',  &
                 ' UNDERLOW CATEGORY CANNOT BE COMPARED.'//  &
                 ' .FALSE. HAS BEEN RETURNED.'/)
         GO TO 170
      ENDIF

!             Check for zero.

      KFLAG = 0
      IF (MA(2).EQ.0) THEN
          JCOMP = 2
          IF (MB(2).LT.0) JCOMP = 1
          IF (MB(2).GT.0) JCOMP = 3
          GO TO 160
      ENDIF
      IF (MB(2).EQ.0) THEN
          JCOMP = 1
          IF (MA(2).LT.0) JCOMP = 3
          GO TO 160
      ENDIF
!             Check for opposite signs.

      IF (MA(2).GT.0 .AND. MB(2).LT.0) THEN
          JCOMP = 1
          GO TO 160
      ENDIF
      IF (MB(2).GT.0 .AND. MA(2).LT.0) THEN
          JCOMP = 3
          GO TO 160
      ENDIF

!             See which one is larger in absolute value.

      IF (MA(1).GT.MB(1)) THEN
          JCOMP = 1
          GO TO 160
      ENDIF
      IF (MB(1).GT.MA(1)) THEN
          JCOMP = 3
          GO TO 160
      ENDIF
      NLAST = NDIG + 1

      DO 150 J = 2, NLAST
         IF (ABS(MA(J)).GT.ABS(MB(J))) THEN
             JCOMP = 1
             GO TO 160
         ENDIF
         IF (ABS(MB(J)).GT.ABS(MA(J))) THEN
             JCOMP = 3
             GO TO 160
         ENDIF
  150 CONTINUE

      JCOMP = 2

!             Now match the JCOMP value to the requested comparison.

  160 IF (JCOMP.EQ.1 .AND. MA(2).LT.0) THEN
          JCOMP = 3
      ELSE IF (JCOMP.EQ.3 .AND. MB(2).LT.0) THEN
          JCOMP = 1
      ENDIF

      FMCOMP = .FALSE.
      IF (JCOMP.EQ.1 .AND. (LREL.EQ.'GT' .OR. LREL.EQ.'GE' .OR.  &
                            LREL.EQ.'NE')) FMCOMP = .TRUE.

      IF (JCOMP.EQ.2 .AND. (LREL.EQ.'EQ' .OR. LREL.EQ.'GE' .OR.  &
                            LREL.EQ.'LE')) FMCOMP = .TRUE.

      IF (JCOMP.EQ.3 .AND. (LREL.EQ.'NE' .OR. LREL.EQ.'LT' .OR.  &
                            LREL.EQ.'LE')) FMCOMP = .TRUE.

  170 IF (NTRACE.NE.0) THEN
          IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.1) THEN
              IF (KFLAG.EQ.0) THEN
                  WRITE (KW,180) NCALL,JBASE,NDIG
  180             FORMAT(' FMCOMP',15X,'CALL LEVEL =',I2,5X,'JBASE =',  &
                         I10,5X,'NDIG =',I6)
              ELSE
                  WRITE (KW,190) NCALL,JBASE,NDIG,KFLAG
  190             FORMAT(' FMCOMP',6X,'CALL LEVEL =',I2,4X,'JBASE =',  &
                         I10,4X,'NDIG =',I6,4X,'KFLAG =',I3)
              ENDIF
              IF (FMCOMP) THEN
                  WRITE (KW,200)
  200             FORMAT(7X,'.TRUE.')
              ELSE
                  WRITE (KW,210)
  210             FORMAT(7X,'.FALSE.')
              ENDIF
          ENDIF
      ENDIF
      NCALL = NCALL - 1
END FUNCTION

SUBROUTINE FMCOS(MA,MB)

!  MB = COS(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(9,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Reduce the argument, convert to radians if the input is
!             in degrees, and evaluate the function.

      CALL FMRDC(MB,MB,JSIN,JCOS,JSWAP)
      IF (MB(1).EQ.KUNKNO) GO TO 110
      IF (KRAD.EQ.0) THEN
          IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMPI2(MPISAV)
              NJBPI = JBASE
              NDIGPI = NDIG
              NDIG = NDSV
          ENDIF
          CALL FMMPY(MB,MPISAV,MB)
          CALL FMDIVI(MB,180,MB)
      ENDIF
      IF (MB(1).NE.KUNKNO) THEN
          IF (JSWAP.EQ.0) THEN
              CALL FMCOS2(MB,MB)
          ELSE
              CALL FMSIN2(MB,MB)
          ENDIF
      ENDIF

!             Append the sign, round, and return.

      IF (JCOS.EQ.-1) MB(2) = -MB(2)
  110 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMCOS2(MA,MB)

!  Internal routine for computing MB = COS(MA) where 0.LE.MA.LE.1.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Use COS(X) = SQRT(1-SIN(X)**2).

      CALL FMSIN2(MA,MB)
      CALL FMMPY(MB,MB,MB)
      CALL FMI2M(1,M02)
      CALL FMSUB(M02,MB,M02)
      CALL FMSQRT(M02,MB)
END SUBROUTINE

SUBROUTINE FMCOSH(MA,MB)

!  MB = COSH(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(10,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      KWRNSV = KWARN
      KWARN = 0
      KARGSW = 0
      CALL FMEQU(MA,MB,NDSAVE,NDIG)
      CALL FMEXP(MB,MB)
      IF (MB(1).EQ.KEXPOV) THEN
          GO TO 110
      ELSE IF (MB(1).EQ.KEXPUN) THEN
          MB(1) = KEXPOV
          GO TO 110
      ENDIF
      CALL FMI2M(1,M01)
      CALL FMDIV(M01,MB,M01)
      CALL FMADD(MB,M01,MB)
      CALL FMDIVI(MB,2,MB)

!             Round and return.

  110 KWARN = KWRNSV
      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMDIG(NSTACK,KST)

!  Compute the number of intermediate digits to be used in Newton
!  iteration.  This assumes that a starting approximation which is
!  accurate to double precision is used, and the root is simple.

!  KST is the number of iterations needed for final accuracy NDIG.
!  NSTACK(J) holds the value of NDIG to be used for the
!            Jth iteration.

      DOUBLE PRECISION XBASE,ONE,Y,YT
      DIMENSION NSTACK(19)

!             Compute the approximate machine precision.  This will be
!             the initial accuracy obtained by using double precision
!             to compute the first approximation before starting Newton
!             iteration.

      XBASE = JBASE
      ONE = 1.0
      Y = ONE
      NE = 1

  110 Y = Y/XBASE
      YT = ONE + Y
      IF (YT.GT.ONE) THEN
          NE = NE + 1
          GO TO 110
      ENDIF

!             Fill the intermediate digit stack (backwards).

      KST = 1
      ND = NDIG
      NSTACK(1) = ND
      IF (ND.LT.2*NE .OR. ND.LE.2) RETURN

  120 Y = ND

!             The 1.9 accounts for the fact that the number of correct
!             digits approximately doubles at each iteration.

      NDT = Y/1.9
      IF (2*NDT.LE.ND) NDT = NDT + 1
      ND = NDT
      KST = KST + 1
      NSTACK(KST) = ND
      IF (ND.GE.2*NE .AND. ND.GT.2) GO TO 120

!             Reverse the stack.

      L = KST/2
      DO 130 J = 1, L
         JT = NSTACK(J)
         NSTACK(J) = NSTACK(KST+1-J)
         NSTACK(KST+1-J) = JT
  130 CONTINUE

END SUBROUTINE

SUBROUTINE FMDIM(MA,MB,MC)

!  MC = DIM(MA,MB)

!  Positive difference.  MC = MA - MB  if MA.GE.MB,
!                           = 0        otherwise.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      CALL FMENTR(11,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN
      KARGSW = 0
      KWSV = KWARN
      KWARN = 0
      MXEXP = MXSAVE

      CALL FMEQU(MA,M01,NDSAVE,NDIG)
      CALL FMEQU(MB,M02,NDSAVE,NDIG)

      IF (FMCOMP(M01,'LT',M02)) THEN
          CALL FMI2M(0,MC)
      ELSE
          CALL FMSUB(M01,M02,MC)
      ENDIF

      IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMDIV(MA,MB,MC)

!  MC = MA / MB

!  This routine performs the trace printing for division.
!  FMDIV2 is used to do the arithmetic.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 12
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      CALL FMDIV2(MA,MB,MC)

      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMDIV2(MA,MB,MC)

!  Internal division routine.  MC = MA / MB

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      IF (KARGSW.NE.1) THEN
          CALL FMARGS(12,2,MA,MB,KRESLT)
          IF (KRESLT.NE.0) THEN
              CALL FMRSLT(MA,MB,MC,KRESLT)
              RETURN
          ENDIF
      ELSE
          IF (MB(2).EQ.0) THEN
              CALL FMIM(0,MC)
              MC(1) = KUNKNO
              MC(2) = 1
              KFLAG = -4
              CALL FMWARN
              RETURN
          ENDIF
          IF (MA(2).EQ.0) THEN
              CALL FMIM(0,MC)
              RETURN
          ENDIF
      ENDIF
      KFLAG = 0

!             NGUARD is the number of guard digits used.

      IF (NCALL.GT.1) THEN
          NGUARD = 1
      ELSE
          B = JBASE
          NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0
      ENDIF
      N1 = NDIG + 1
      NG = N1 + NGUARD

!             Copy MA into the working array.

      DO 110 J = 3, N1
         MWA(J+1) = MA(J)
  110 CONTINUE
      MWA(1) = MA(1) - MB(1) + 1
      MWA(2) = 0
      N3 = NDIG + 3
      NL = N1 + NGUARD + 2
      DO 120 J = N3, NL
         MWA(J) = 0
  120 CONTINUE

!             Save the sign of MA and MB and then work only with
!             positive numbers.

      MA2 = MA(2)
      MB1 = MB(1)
      MB2 = MB(2)
      MA(2) = ABS(MA(2))
      MWA(3) = MA(2)
      MB(1) = 0
      MB(2) = ABS(MB(2))

!             NMBWDS is the number of words of MB used to compute
!             the trial quotient digit LT.

      NMBWDS = 2
      IF (JBASE.LT.100) NMBWDS = 4

!             XB is an approximation of MB used in selecting
!             trial quotients.

      XBASE = JBASE
      XB = 0
      JL = NMBWDS + 1
      DO 130 J = 2, JL
         IF (J.LE.N1) THEN
             XB = XB*XBASE + MB(J)
         ELSE
             XB = XB*XBASE
         ENDIF
  130 CONTINUE
      IF (JL+1.LE.N1) XB = XB + MB(JL+1)/XBASE
      XBLOG = LOG(XBASE)

!             LMAX determines when to normalize all of MWA.
!             NGAP determines how many digits must be normalized for
!             each trial quotient.

      IF (MXBASE**2/(JBASE-1)**2.GE.NL) THEN
          LMAX = NL*(JBASE-1)
      ELSE
          LMAX = MXBASE**2/(JBASE-1)
      ENDIF
      X = LMAX
      NGAP = LOG(X)/XBLOG + 2.0

!             Count the trailing zero digits of MB.

      DO 140 J = 1, NDIG
         IF (MB(NDIG+2-J).NE.0) THEN
             NZDMB = J - 1
             GO TO 150
         ENDIF
  140 CONTINUE

!             MAXMWA is an upper bound on the size of values in MWA
!             divided by JBASE-1.  It is used to determine whether
!             normalization can be postponed.

  150 MAXMWA = 0

!             KPTMWA points to the next digit in the quotient.

      KPTMWA = 2

!             This is the start of the division loop.

!             LT is the trial quotient digit.
!             Divide the first (NMBWDS+1) digits of MWA by MB so that
!             LT has a high probability of being the correct value.

  160 KL = KPTMWA + NMBWDS
      XMWA = 0
      DO 170 J = KPTMWA, KL
         IF (J.LE.NL) THEN
             XMWA = XMWA*XBASE + MWA(J)
         ELSE
             XMWA = XMWA*XBASE
         ENDIF
  170 CONTINUE

      LT = XMWA/XB
      IF (LT.GE.JBASE) LT = JBASE - 1

!             Subtract LT*MB from MWA.

      KA = KPTMWA + 1
      KB = MIN(KA+NDIG-1-NZDMB,NL)
      JB = KA - 2
      IF (LT.GT.0) THEN
          DO 180 J = KA, KB
             MWA(J) = MWA(J) - LT*MB(J-JB)
  180     CONTINUE
      ENDIF

!             Normalize enough digits in MWA so that LQ, the correct
!             quotient digit, can be determined.

      LQ = LT
      MAXMWA = MAXMWA + LT

!             NORMPT points to the last digit currently normalized.

      NORMPT = MIN(KA+NGAP+1,KB)

!             If the correct LQ has been found, then the current section
!             of MWA will be nonnegative and less than MB.
!             If MWA.LT.0 then too large a trial quotient was used and
!             a '+ MB' correction step is needed.
!             If MWA.GE.MB then too small a trial quotient was used and
!             a '- MB' correction step is needed.

  190 J = NORMPT
  200 IF (MWA(J).LT.0) THEN
          KT = (-MWA(J)-1)/JBASE + 1
          MWA(J-1) = MWA(J-1) - KT
          MWA(J) = MWA(J) + KT*JBASE
      ENDIF
      J = J - 1
      IF (J.GE.KA) GO TO 200
      IF (NORMPT.EQ.KB) MAXMWA = 0

!             NONZPT points to the first nonzero digit in MWA after the
!             subtraction and normalization have been done.

      IF (MWA(KA-1).NE.0) THEN
          NONZPT = KA - 1
      ELSE IF (MWA(KA).NE.0) THEN
          NONZPT = KA
      ELSE
          NONZPT = 0
          K1 = KA + 1
          IF (K1.LE.NORMPT) THEN
              DO 210 J = K1, NORMPT
                 IF (MWA(J).NE.0) THEN
                     NONZPT = J
                     GO TO 220
                 ENDIF
  210         CONTINUE
          ENDIF
      ENDIF

!             KDIFF points to the first digit in MWA which differs
!             from the corresponding digit in MB.

  220 IF (MWA(KA-1).NE.0) THEN
          KDIFF = KA - 1
      ELSE IF (MWA(KA).NE.MB(2)) THEN
          KDIFF = KA
      ELSE
          KDIFF = 0
          K1 = KA + 1
          IF (K1.LE.NORMPT) THEN
              DO 230 J = K1, NORMPT
                 IF (MWA(J).NE.MB(J-JB)) THEN
                     KDIFF = J
                     GO TO 240
                 ENDIF
  230         CONTINUE
          ENDIF
      ENDIF

!             Protect against cases like MWA = 2 0 0 0 ...
!             and MB = 1 9 9 9 ....  Since MWA may not be completely
!             normalized, MWA(KDIFF) might still be equal to
!             MB(KDIFF-JB).

  240 IF (KDIFF.NE.0 .AND. KDIFF.LT.NL .AND. KDIFF+1-JB.LE.N1) THEN
          IF (MWA(KDIFF+1).EQ.0 .AND. NORMPT.LT.KB .AND.  &
              MB(KDIFF+1-JB).EQ.JBASE-1) KDIFF = 0
      ENDIF

!             If NONZPT is zero then not enough digits were normalized
!             to tell whether MWA is negative.

      IF (NONZPT.EQ.0 .OR. NORMPT-NONZPT.LT.NGAP) THEN
          IF (NORMPT.EQ.KB) GO TO 250
          NORMPT = KB
          GO TO 190
      ENDIF

!             See if the current section of MWA is negative.

  250 IF (NONZPT.LE.0) GO TO 270
      IF (MWA(NONZPT).LT.0) THEN

!             At the end of the division LT could be off by more than
!             one because part of MB has shifted off the end of the
!             active part of MWA.  Skip the correction step for the
!             last guard digit.

          IF (KA.EQ.KB .AND. KB.EQ.NL) GO TO 310

!             Correct LQ by adding MB.

          LQ = LQ - 1
          MAXMWA = MAXMWA - 1

!             This is the only place where values greater than JBASE
!             can be introduced into MWA.  Normalize them here.

          K = KB
          KC = 0
  260     MWA(K) = MWA(K) + MB(K-JB) + KC
          IF (MWA(K).GE.JBASE) THEN
              KC = 1
              MWA(K) = MWA(K) - JBASE
          ELSE
              KC = 0
          ENDIF
          K = K - 1
          IF (K.GE.KA) GO TO 260
          IF (KC.EQ.1) MWA(KA-1) = MWA(KA-1) + 1
          GO TO 190
      ENDIF

!             If KDIFF is zero then not enough digits were normalized
!             to tell if MWA.GE.MB.

  270 IF (KDIFF.EQ.0) THEN
          IF (NORMPT.EQ.KB) GO TO 280
          NORMPT = KB
          GO TO 190
      ENDIF

!             See if the current section of MWA is greater than MB.

      IF (MWA(KDIFF).LT.MB(KDIFF-JB)) GO TO 310

!             If NORMPT-KDIFF.GE.NGAP then not enough digits were
!             normalized to tell if MWA.GE.MB.

      IF (NORMPT-KDIFF.LT.NGAP .AND. NORMPT.LT.KB) THEN
          NORMPT = KB
          GO TO 190
      ENDIF

!             Correct LQ by subtracting MB.

  280 IF (KA.EQ.KB .AND. KB.EQ.NL) GO TO 310
      LQ = LQ + 1

!             Due to part of MB shifting off the end of the active
!             part of MWA, LQ could be equal to JBASE here.  This
!             means the correct value of LQ is JBASE-1 for all
!             subsequent digits.

      IF (LQ.GE.JBASE) THEN
          NG1 = NG + 1
          DO 290 K = KPTMWA, NG1
             MWA(K) = JBASE - 1
  290     CONTINUE
          GO TO 330
      ENDIF

      MAXMWA = MAXMWA + 1
      DO 300 K = KA, KB
         MWA(K) = MWA(K) - MB(K-JB)
  300 CONTINUE
      GO TO 190

!             Here the correct value for LQ has been found.

  310 MWA(KPTMWA) = LQ

!             See if MWA must be normalized.

      IF (MAXMWA+JBASE-1.GE.LMAX) THEN
          J = KB
  320     IF (MWA(J).LT.0) THEN
              KT = (-MWA(J)-1)/JBASE + 1
              MWA(J-1) = MWA(J-1) - KT
              MWA(J) = MWA(J) + KT*JBASE
          ENDIF
          J = J - 1
          IF (J.GE.KA) GO TO 320
          MAXMWA = 0
      ENDIF

      KPTMWA = KPTMWA + 1
      IF (KPTMWA.LE.NG) GO TO 160
      IF (MWA(2).EQ.0 .AND. KPTMWA.LE.NG+1) GO TO 160

!             Round, affix the sign, and return.

  330 MA(2) = MA2
      MB(1) = MB1
      MB(2) = MB2
      IF (MWA(2).EQ.0) THEN
          CALL FMRND(NDIG,NGUARD,1)
      ELSE
          CALL FMRND(NDIG,NGUARD,0)
      ENDIF
      CALL FMMOVE(MC)
      IF (MA2*MB2.LT.0) MC(2) = -MC(2)
END SUBROUTINE

SUBROUTINE FMDIVI(MA,INT_T,MB)

!  MB = MA / INT_T

!  Divide FM number MA by one word integer INT_T.

!  This routine is faster than FMDIV when the divisor is less than
!  MXBASE.  When INT_T is not less than MXBASE, FMDIV2 is used.  In
!  this case if INT_T is known to be a product of two integers less than
!  MXBASE it is usually faster to make two calls to FMDIVI with half-
!  word factors than one call with their product.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 13
      IF (NTRACE.NE.0) THEN
          CALL FMNTR(2,MA,MA,1)
          CALL FMNTRI(2,INT_T,0)
      ENDIF
      KFLAG = 0
      N1 = NDIG + 1

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO .OR. INT_T.EQ.0) THEN
          MA1 = MA(1)
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          IF (MA1.NE.KUNKNO) CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (ABS(INT_T).EQ.1) THEN
          DO 110 J = 1, N1
             MB(J) = MA(J)
  110     CONTINUE
          MB(2) = MA(2)*INT_T
          IF (MA(1).EQ.KEXPOV) KFLAG = -5
          IF (MA(1).EQ.KEXPUN) KFLAG = -6
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(2).EQ.0) THEN
          DO 120 J = 1, N1
             MB(J) = 0
  120     CONTINUE
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPUN) THEN
          MA2 = MA(2)
          CALL FMIM(0,MB)
          MB(1) = KEXPUN
          MB(2) = 1
          IF ((MA2.LT.0 .AND. INT_T.GT.0) .OR.  &
              (MA2.GT.0 .AND. INT_T.LT.0))      MB(2) = -1
          KFLAG = -6
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPOV) THEN
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             NGUARD is the number of guard digits used.

      IF (NCALL.GT.1) THEN
          NGUARD = 1
      ELSE
          B = JBASE
          NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0
      ENDIF
      N1 = NDIG + 1

!             If ABS(INT_T).GE.MXBASE use FMDIV.

      IF (ABS(INT_T).GT.MXBASE) THEN
          CALL FMIM(INT_T,M01)
          CALL FMDIV2(MA,M01,MB)
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Work with positive numbers.

      MA2 = MA(2)
      MA(2) = ABS(MA(2))
      INTP = ABS(INT_T)

!             Find the first significant digit of the quotient.

      KT = 0
      DO 130 J = 2, N1
         KPT = J
         KT = KT*JBASE + MA(J)
         IF (KT.GE.INTP) GO TO 150
  130 CONTINUE

  140 KPT = KPT + 1
      KT = KT*JBASE
      IF (KT.LT.INTP) GO TO 140

!             Do the rest of the division.

  150 KA = KPT + 1
      MWA(1) = MA(1) + 2 - KPT
      MWA(2) = KT/INTP
      MODINT = MOD(KT,INTP)
      KPTWA = 2
      IF (KA.LE.N1) THEN
          KL = 3 - KA
          DO 160 J = KA, N1
             KT = MODINT*JBASE + MA(J)
             MWA(KL+J) = KT/INTP
             MODINT = MOD(KT,INTP)
  160     CONTINUE
          KPTWA = KL + N1
      ENDIF

      KA = KPTWA + 1
      KB = N1 + NGUARD
      IF (KA.LE.KB) THEN
          DO 170 J = KA, KB
             KT = MODINT*JBASE
             MWA(J) = KT/INTP
             MODINT = MOD(KT,INTP)
  170     CONTINUE
      ENDIF

!             Round the result, put the sign on MB and return.

      MA(2) = MA2
      CALL FMRND(NDIG,NGUARD,0)
      CALL FMMOVE(MB)
      IF ((MA2.LT.0 .AND. INT_T.GT.0) .OR. (MA2.GT.0 .AND. INT_T.LT.0))  &
           MB(2) = -MB(2)
      IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMDM(X,MA)

!  Internal routine for converting double precision to multiple
!  precision.

      DOUBLE PRECISION ONE,X,XBASE,Y,YT
      DIMENSION MA(LUNPCK)

      KFLAG = 0
      N1 = NDIG + 1
      ONE = 1.0
      XBASE = JBASE
      Y = ONE
      K = 0
      NE = 2

!             Find the approximate machine precision.
!             NE-1 is the number of words at the current precision and
!             base roughly equal to machine precision.

      DO 110 J = 2, NDIG
         Y = Y/XBASE
         NE = NE + 1
         YT = ONE + Y
         IF (YT.LE.ONE) GO TO 120
  110 CONTINUE
  120 Y = X
      IF (X.LT.0.0) Y = -X

      IF (X.EQ.0.0) THEN
          DO 130 J = 1, N1
             MA(J) = 0
  130     CONTINUE
          GO TO 240
      ENDIF

!             Get the exponent.

      IF (Y.GT.ONE) THEN
          IF (Y/XBASE.LT.Y) THEN
  140         K = K + 1
              Y = Y/XBASE
              IF (Y.GT.ONE) GO TO 140
              IF (Y.LT.ONE) THEN
                  MA(1) = K
                  GO TO 180
              ENDIF
              GO TO 160
          ELSE
              CALL FMIM(0,MA)
              MA(1) = KUNKNO
              MA(2) = 1
              KFLAG = -4
              CALL FMWARN
              RETURN
          ENDIF
      ENDIF

      IF (Y.LT.ONE) THEN
          IF (Y*XBASE.GT.Y) THEN
  150         K = K - 1
              Y = Y*XBASE
              IF (Y.LT.ONE) GO TO 150
              IF (Y.GT.ONE) THEN
                  K = K + 1
                  Y = Y/XBASE
                  MA(1) = K
                  GO TO 180
              ENDIF
          ELSE
              CALL FMIM(0,MA)
              MA(1) = KUNKNO
              MA(2) = 1
              KFLAG = -4
              CALL FMWARN
              RETURN
          ENDIF
      ENDIF

  160 MA(1) = K + 1
      MA(2) = 1
      DO 170 J = 3, N1
         MA(J) = 0
  170 CONTINUE
      GO TO 240

!             Build the rest of the number.

  180 DO 190 J = 2, NE
         Y = Y*XBASE
         K = Y
         YT = K
         Y = Y - YT
         MA(J) = K
         IF (J.GE.N1) GO TO 210
  190 CONTINUE
      K = NE + 1
      DO 200 J = K, N1
         MA(J) = 0
  200 CONTINUE

!             Normalize.

  210 IF (ABS(MA(2)).GE.JBASE) THEN
          K = N1 + 1
          DO 220 J = 3, N1
             K = K - 1
             MA(K) = MA(K-1)
  220     CONTINUE
          MN = MA(2)/JBASE
          MA(3) = MA(2) - MN*JBASE
          MA(2) = MN
          MA(1) = MA(1) + 1
          GO TO 240
      ENDIF

      IF (MA(2).EQ.0) THEN
          DO 230 J = 2, NDIG
             MA(J) = MA(J+1)
  230     CONTINUE
          MA(1) = MA(1) - 1
          MA(N1) = 0
      ENDIF

  240 IF (X.LT.0.0) MA(2) = -MA(2)
END SUBROUTINE

SUBROUTINE FMDP2M(X,MA)

!  MA = X

!  Convert a double precision floating point number to FM format.

!  In general the relative accuracy of the number returned is only
!  the relative accuracy of a machine precision number.  This may be
!  true even if X can be represented exactly in the machine floating
!  point number system.

      DOUBLE PRECISION X
      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 14
      IF (NTRACE.NE.0) CALL FMNTRR(2,X,1)

      CALL FMDM(X,MA)

      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMENTR(NROUTN,MA,MB,NARGS,MC,KRESLT,NDSAVE,MXSAVE,  &
                        KASAVE,KOVUN)

!  Do the argument checking and increasing of precision, overflow
!  threshold, etc., upon entry to an FM routine.

!  NROUTN - routine number of calling routine
!  MA     - first input argument
!  MB     - second input argument (optional)
!  NARGS  - number of input arguments
!  MC     - result argument
!  KRESLT - returned nonzero if the input arguments give the result
!           immediately (e.g., MA*0 or OVERFLOW*MB)
!  NDSAVE - saves the value of NDIG after NDIG is increased
!  MXSAVE - saves the value of MXEXP
!  KASAVE - saves the value of KARGSW
!  KOVUN  - returned nonzero if an input argument is (+ or -) overflow
!           or underflow.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = NROUTN
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,NARGS)
      CALL FMARGS(NROUTN,NARGS,MA,MB,KRESLT)

      KOVUN = 0
      IF (MA(1).EQ.KEXPOV .OR. MA(1).EQ.KEXPUN) KOVUN = 1
      IF (NARGS.EQ.2) THEN
          IF (MB(1).EQ.KEXPOV .OR. MB(1).EQ.KEXPUN) KOVUN = 1
      ENDIF

!             Increase the working precision.

      NDSAVE = NDIG
      IF (NCALL.EQ.1) THEN
          B = JBASE
          K = 5.0*LOG(10.0)/LOG(B) + 2.0
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              KRESLT = 12
          ENDIF
      ENDIF

      IF (KRESLT.NE.0) THEN
          IF (KRESLT.EQ.9 .OR. KRESLT.EQ.10 .OR. KRESLT.GE.13) THEN
              IF (KRAD.EQ.1) THEN
                  CALL FMPI(MC)
              ELSE
                  CALL FMI2M(180,MC)
              ENDIF
              IF (KRESLT.LE.10) CALL FMDIVI(MC,2,MC)
              IF (KRESLT.GE.14) CALL FMDIVI(MC,4,MC)
              CALL FMEQU(MC,MC,NDIG,NDSAVE)
              NDIG = NDSAVE
              IF (KRESLT.EQ.9 .OR. KRESLT.EQ.14) MC(2) = -MC(2)
              IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
              NCALL = NCALL - 1
              RETURN
          ENDIF

          NDIG = NDSAVE
          CALL FMRSLT(MA,MB,MC,KRESLT)
          IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Turn off argument checking for low-level routines while
!             a high-level call is in progress.

      KASAVE = KARGSW
      KARGSW = 1

!             Extend the overflow/underflow threshold.

      MXSAVE = MXEXP
      MXEXP = MXEXP2
END SUBROUTINE

SUBROUTINE FMEQU(MA,MB,NDA,NDB)

!  Set MB (having NDB digits) equal to MA (having NDA digits).

!  If MB has less precision than MA the result is rounded to NDB digits.

!  If MB has more precision the result has zero digits padded on the
!  right.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             See if MA and MB are the same array.  If so, the loop
!             which copies MA to MB can be skipped.

!             The MIN calls try to obscure the test so that an
!             optimizing compiler will not remove the K = ... and
!             MA(3) = ... lines and sometimes change an input argument.

      KSAME = 0
      K = MIN(MA(3),MXBASE)
      MB(3) = K + 1
      IF (MA(3).EQ.MB(3)) KSAME = 1
      MA(3) = MIN(K,MXBASE+MXEXP)

!             Check for special symbols.

      KFLAG = 0
      IF (ABS(MA(1)).GE.KEXPOV) THEN
          DO 110 J = 2, NDB
             MB(J+1) = 0
  110     CONTINUE
          MB(1) = MA(1)
          MB(2) = MA(2)
          GO TO 240
      ENDIF

      IF (NDB.EQ.NDA) GO TO 190

      IF (NDB.GT.NDA) GO TO 210

!             Round to NDB digits.

      NDG = NDB
      N1 = NDB + 1
      IF (KSAME.EQ.0) THEN
          DO 120 J = 1, N1
             MB(J) = MA(J)
  120     CONTINUE
      ENDIF
      IF (NDG.LT.1 .OR. (KROUND.EQ.0 .AND. NCALL.LE.1)) GO TO 240

      L = NDB + 2
      IF (2*(MA(L)+1).LT.JBASE) GO TO 240
      IF (MOD(JBASE,2).EQ.0) THEN
          IF (2*MA(L).LT.JBASE) GO TO 240
          IF (2*MA(L).EQ.JBASE) THEN
              IF (L.LE.NDA) THEN
                  DO 130 J = L, NDA
                     IF (MA(J+1).GT.0) GO TO 150
  130             CONTINUE
              ENDIF

!                       Round to even.

              IF (MOD(MB(N1),2).EQ.0) GO TO 240
          ENDIF
      ELSE
          IF (2*MA(L)+1.EQ.JBASE) THEN
              IF (L.LE.NDA) THEN
                  DO 140 J = L, NDA
                     IF (2*(MA(J+1)+1).LT.JBASE) GO TO 240
                     IF (2*MA(J+1).GT.JBASE) GO TO 150
  140             CONTINUE
                  GO TO 240
              ENDIF
          ENDIF
      ENDIF

  150 MB(NDG+1) = MB(NDG+1) + 1
      MB(NDG+2) = 0

!             Check whether there was a carry in the rounded digit.

      MB2 = MB(2)
      MB(2) = ABS(MB(2))
      KB = NDG + 1
      IF (KB.GE.3) THEN
          K = KB + 1
          DO 160 J = 3, KB
             K = K - 1
             IF (MB(K).LT.JBASE) GO TO 180
             MB(K-1) = MB(K-1) + MB(K)/JBASE
             MB(K) = MOD(MB(K),JBASE)
  160     CONTINUE
      ENDIF

!             If there is a carry in the first digit then the exponent
!             must be adjusted and the number shifted right.

      IF (MB(2).LT.JBASE) GO TO 180
      IF (KB.GE.4) THEN
          K = KB + 1
          DO 170 J = 4, KB
             K = K - 1
             MB(K) = MB(K-1)
  170     CONTINUE
      ENDIF

      IF (KB.GE.3) MB(3) = MOD(MB(2),JBASE)
      MB(2) = MB(2)/JBASE
      MB(1) = MB(1) + 1

  180 IF (MB2.LT.0) MB(2) = -MB(2)
      GO TO 240

!             MA and MB have the same precision.

  190 IF (KSAME.EQ.0) THEN
          NDA1 = NDA + 1
          DO 200 J = 1, NDA1
             MB(J) = MA(J)
  200     CONTINUE
      ENDIF
      GO TO 240

!             Extend to NDB digits by padding with zeros.

  210 IF (KSAME.EQ.0) THEN
          NDA1 = NDA + 1
          DO 220 J = 1, NDA1
             MB(J) = MA(J)
  220     CONTINUE
      ENDIF
      NDA2 = NDA + 2
      NDB2 = NDB + 1
      DO 230 J = NDA2, NDB2
         MB(J) = 0
  230 CONTINUE

!             Check for overflow or underflow.

  240 IF (ABS(MB(1)).GT.MXEXP) THEN
          IF (MB(1).NE.KUNKNO .OR. MB(2).NE.1) THEN
              KWRNSV = KWARN
              KWARN = 0
              NCALL = NCALL + 1
              CALL FMTRAP(MB)
              NCALL = NCALL - 1
              KWARN = KWRNSV
          ENDIF
          IF (MB(1).EQ.KUNKNO) KFLAG = -4
      ENDIF

END SUBROUTINE

SUBROUTINE FMEXIT(MT,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)

!  Upon exit from an FM routine the result MT (having precision NDIG)
!  is rounded and returned in MC (having precision NDSAVE).
!  The values of NDIG, MXEXP, and KARGSW are restored.
!  KOVUN is nonzero if one of the routine's input arguments was overflow
!  or underflow.

      DIMENSION MT(LUNPCK),MC(LUNPCK)

      KWSV = KWARN
      KWARN = 0
      MXEXP = MXSAVE
      KFSAVE = KFLAG
      CALL FMEQU(MT,MC,NDIG,NDSAVE)
      IF (KFLAG.NE.-5 .AND. KFLAG.NE.-6) KFLAG = KFSAVE
      NDIG = NDSAVE
      KWARN = KWSV
      IF (KFLAG.EQ.1) KFLAG = 0
      IF ((MC(1).EQ.KUNKNO .AND. KFLAG.NE.-9)  &
         .OR. (MC(1).EQ.KEXPUN .AND. KOVUN.EQ.0)  &
         .OR. (MC(1).EQ.KEXPOV .AND. KOVUN.EQ.0)) CALL FMWARN
      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
      KARGSW = KASAVE
END SUBROUTINE

SUBROUTINE FMEXP(MA,MB)

!  MB = EXP(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(15,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      MA1 = MA(1)
      MA2 = MA(2)
      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Check for obvious underflow or overflow.
!             XOV is LN(LN(slightly above overflow))
!             XMA is LN(LN(EXP(MA))) approximately.

      XOV = LOG(1.01*REAL(MXEXP)) + LOG(LOG(REAL(JBASE)))
      XMA = LOG(REAL(MAX(ABS(MA2),1))) - LOG(REAL(JBASE)) +  &
            MA1*LOG(REAL(JBASE))

  110 IF (XMA.GE.XOV) THEN
          CALL FMIM(0,MB)
          IF (MA2.GT.0) THEN
              KFLAG = -5
              MB(1) = KEXPOV
              MB(2) = 1
          ELSE
              KFLAG = -6
              MB(1) = KEXPUN
              MB(2) = 1
          ENDIF
          NDIG = NDSAVE
          MXEXP = MXSAVE
          KARGSW = KASAVE
          CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Split MA into integer and fraction parts.
!             Work with a positive argument.
!             M02 = integer part of ABS(MA)
!             MB  = fraction part of ABS(MA)

      MB(2) = ABS(MB(2))
      CALL FMINT(MB,M02)
      CALL FMSUB(MB,M02,MB)

!             If the integer part is not zero, use FMIPWR to compute
!             E**(M02).  If M02 is too large to represent as a one word
!             integer, the definition of MXEXP insures that E**(M02)
!             overflows or underflows.

      KWSV = KWARN
      KWARN = 0
      CALL FMM2I(M02,KT)
      KWARN = KWSV
      IF (KFLAG.NE.0) THEN
          XMA = XOV
          GO TO 110
      ENDIF
      IF (KT.GT.0) THEN

!             Compute IEXTRA, the number of extra digits required in
!             order to get EXP(KT) correct to the current precision.
!             IEXTRA = (LN(KT) - 1)/LN(JBASE)

          XT = KT
          XBASE = JBASE
          IEXTRA = LOG(XT)/LOG(XBASE) + 0.5
          IF (IEXTRA.GT.0 .AND. NDIG+IEXTRA.LE.MXNDG2) THEN
              DO 120 J = 1, IEXTRA
                 MB(NDIG+1+J) = 0
  120         CONTINUE
          ENDIF
          NDIG = NDIG + IEXTRA
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MB(1) = KUNKNO
              MB(2) = 1
              DO 130 J = 2, NDSAVE
                 MB(J+1) = 0
  130         CONTINUE
              CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
              RETURN
          ENDIF

!             Check whether the current precision of e is large
!             enough.

          IF (NJBE.NE.JBASE .OR. NDIG.GT.NDIGE) THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMI2M(1,MESAV)
              CALL FMEXP2(MESAV,MESAV)
              NJBE = JBASE
              NDIGE = NDIG
              NDIG = NDSV
          ENDIF

      ENDIF

!             Now do the fraction part of MA and combine the results.

      KWSV = KWARN
      KWARN = 0
      IF (MB(2).NE.0 .AND. KT.GT.0) THEN
          CALL FMEXP2(MB,MB)
          CALL FMIPWR(MESAV,KT,M03)
          CALL FMMPY(MB,M03,MB)
      ELSE IF (MB(2).NE.0 .AND. KT.EQ.0) THEN
          CALL FMEXP2(MB,MB)
      ELSE IF (MB(2).EQ.0 .AND. KT.GT.0) THEN
          CALL FMIPWR(MESAV,KT,MB)
      ELSE
          CALL FMI2M(1,MB)
      ENDIF

!             Invert if MA was negative.

      IF (MA2.LT.0) THEN
          CALL FMI2M(1,M02)
          CALL FMDIV(M02,MB,MB)
      ENDIF
      KWARN = KWSV

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMEXP2(MA,MB)

!  MB = EXP(MA)

!  Internal exponential routine (called with 0.LT.MA.LE.1).

!  Upon return M03 (in COMMON /FM1/) contains an accurate value for
!  EXP(MA)-1 which can be used to avoid loss of significance in other
!  routines if MA is close to zero.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             LJSUMS = 8*LUNPCK allows for up to eight concurrent sums.
!             Increasing this value will begin to improve the speed of
!             EXP when the base is large and precision exceeds about
!             1,500 decimal digits.

      NDSAVE = NDIG
      IF (MA(1).EQ.1) THEN

!             Here the special case EXP(1.0) is computed.
!             Use the direct series  e = 1/0! + 1/1! + 1/2! + ...
!             Do as much of the work as possible using small integers
!             to minimize the number of FM calls.
!             Reduce NDIG while computing each term in the
!             sum as the terms get smaller.

          B = JBASE
          T = NDIG
          XN = T*LOG(B)/LOG(T)
          K = LOG(XN)/LOG(B)
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MB(1) = KUNKNO
              MB(2) = 1
              DO 110 J = 2, NDSAVE
                 MB(J+1) = 0
  110         CONTINUE
              NDIG = NDSAVE
              RETURN
          ENDIF
          NDSAV1 = NDIG

          CALL FMI2M(2,MB)
          CALL FMI2M(1,M02)
          J = 2
          NBIG = MXBASE

  120     NTOP = 1
          NBOT = J
  130     IF (NBOT.GT.NBIG/(J+1)) GO TO 140
          J = J + 1
          NTOP = J*NTOP + 1
          NBOT = J*NBOT
          GO TO 130

  140     CALL FMDIVI(M02,NBOT,M02)
          IF (NTOP.GT.1) THEN
              CALL FMMPYI(M02,NTOP,M03)
              NDIG = NDSAV1
              CALL FMADD(MB,M03,MB)
              NDIG = NDSAV1 - (MB(1)-M03(1))
          ELSE
              NDIG = NDSAV1
              CALL FMADD(MB,M02,MB)
              NDIG = NDSAV1 - (MB(1)-M02(1))
          ENDIF
          IF (NDIG.LT.2) NDIG = 2
          IF (KFLAG.NE.1) THEN
              J = J + 1
              GO TO 120
          ENDIF
          NDIG = NDSAVE
          CALL FMI2M(-1,M02)
          CALL FMADD(MB,M02,M03)
          KFLAG = 0
          RETURN
      ENDIF

!             Here is the general case.  Compute EXP(MA) where
!             0 .LT. MA .LT. 1.

!             Use the direct series
!                  EXP(X) = 1 + X + X**2/2! + X**3/3! + ...

!             The argument will be halved K2 times before the series
!             is summed.  The series will be added as J2 concurrent
!             series.  The approximately optimal values of K2 and J2
!             are now computed to try to minimize the time required.
!             N2 is the approximate number of terms of the series which
!             will be needed, and L2 guard digits will be carried.

      B = JBASE
      K = 5.0*LOG(10.0)/LOG(B) + 2.0
      T = MAX(NDIG-K,2)
      ALOGB = LOG(B)
      ALOG2 = LOG(2.0)
      ALOGT = LOG(T)
      J2 = 0.081*ALOGB*T**0.3333 + 1.85
      J2 = MAX(1,MIN(J2,LJSUMS/MXNDG2))
      K2 = 0.7*SQRT(T*ALOGB/REAL(J2)) - 0.57*ALOGT + 0.5

      L = -(REAL(MA(1))*ALOGB+LOG(REAL(MA(2))/B +  &
            REAL(MA(3))/(B*B)))/ALOG2 - 0.3
      K2 = K2 - L
      IF (K2.LT.0) THEN
          K2 = 0
          J2 = .43*SQRT(T*ALOGB/(ALOGT+REAL(L)*ALOG2)) + .33
      ENDIF
      IF (J2.LE.1) J2 = 1

      N2 = T*ALOGB/(ALOGT+REAL(L)*ALOG2)
      L2 = LOG(REAL(N2)+2.0**K2)/ALOGB
      NDIG = NDIG + L2
      IF (NDIG.GT.MXNDG2) THEN
          KFLAG = -9
          CALL FMWARN
          MB(1) = KUNKNO
          MB(2) = 1
          DO 150 J = 2, NDSAVE
             MB(J+1) = 0
  150     CONTINUE
          NDIG = NDSAVE
          RETURN
      ENDIF
      NDSAV1 = NDIG

!             Halve the argument K2 times.

      CALL FMEQU(MA,M02,NDSAVE,NDIG)
      KTWO = 1
      MAXVAL = MXBASE/2
      IF (K2.GT.0) THEN
          DO 160 J = 1, K2
             KTWO = 2*KTWO
             IF (KTWO.GT.MAXVAL) THEN
                 CALL FMDIVI(M02,KTWO,M02)
                 KTWO = 1
             ENDIF
  160     CONTINUE
          IF (KTWO.GT.1) CALL FMDIVI(M02,KTWO,M02)
      ENDIF

!             Sum the series X + X**2/2! + X**3/3! + ....
!             Split into J2 concurrent sums and reduce NDIG while
!             computing each term in the sum as the terms get smaller.

      CALL FMEQU(M02,MB,NDIG,NDIG)
      NTERM = 1
      DO 170 J = 1, J2
         CALL FMDIVI(MB,NTERM,MB)
         NTERM = NTERM + 1
         KPT = (J-1)*(NDIG+1) + 1
         CALL FMEQU(MB,MJSUMS(KPT),NDIG,NDIG)
  170 CONTINUE
      CALL FMIPWR(M02,J2,M03)

  180 CALL FMMPY(MB,M03,MB)
      DO 190 J = 1, J2
         CALL FMDIVI(MB,NTERM,MB)
         KPT = (J-1)*(NDSAV1+1) + 1
         NDIG = NDSAV1
         CALL FMADD(MJSUMS(KPT),MB,MJSUMS(KPT))
         IF (KFLAG.EQ.1) GO TO 200
         NDIG = NDSAV1 - (MJSUMS(KPT)-MB(1))
         IF (NDIG.LT.2) NDIG = 2
         NTERM = NTERM + 1
  190 CONTINUE
      GO TO 180

!             Next put the J2 separate sums back together.

  200 KFLAG = 0
      KPT = (J2-1)*(NDIG+1) + 1
      CALL FMEQU(MJSUMS(KPT),M03,NDIG,NDIG)
      IF (J2.GE.2) THEN
          DO 210 J = 2, J2
             CALL FMMPY(M02,M03,M03)
             KPT = (J2-J)*(NDIG+1) + 1
             CALL FMADD(M03,MJSUMS(KPT),M03)
  210     CONTINUE
      ENDIF

!             Reverse the effect of halving the argument to
!             compute EXP(MA).

      NDIG = NDSAV1
      IF (K2.GT.0) THEN
          CALL FMI2M(2,MB)
          DO 220 J = 1, K2
             CALL FMADD(M03,MB,M02)
             CALL FMMPY(M03,M02,M03)
  220     CONTINUE
      ENDIF

      CALL FMI2M(1,MB)
      CALL FMADD(M03,MB,MB)
      CALL FMEQU(MB,MB,NDSAV1,NDSAVE)
      NDIG = NDSAVE

END SUBROUTINE

SUBROUTINE FMI2M(INTEG,MA)

!  MA = INTEG

!  Convert an integer to FM format.

!  The conversion is exact if INTEG is less than JBASE**NDIG,
!  otherwise the result is an approximation.

!  This routine performs the trace printing for the conversion.
!  FMIM is used to do the arithmetic.

      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 16
      IF (NTRACE.NE.0) CALL FMNTRI(2,INTEG,1)

      CALL FMIM(INTEG,MA)

      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMIM(INTEG,MA)

!  MA = INTEG.  Internal integer conversion routine.

!  The conversion is exact if INTEG is less than JBASE**NDIG.
!  Otherwise FMDM is used to get an approximation.

      DIMENSION MA(LUNPCK)
      DOUBLE PRECISION X

      KFLAG = 0
      N1 = NDIG + 1

!              Check for INTEG equal to zero.

      INT_T = ABS(INTEG)
      IF (INT_T.EQ.0) THEN
          DO 110 J = 1, N1
             MA(J) = 0
  110     CONTINUE
          GO TO 150
      ENDIF

!             Compute and store the digits, right to left.

      MA(1) = 0
      J = NDIG + 1

  120 K = INT_T/JBASE
      L = INT_T - K*JBASE
      MA(1) = MA(1) + 1
      MA(J) = L
      IF (K.GT.0) THEN
          INT_T = K
          J = J - 1
          IF (J.GE.2) GO TO 120

!             Here INTEG cannot be expressed exactly.

          X = INTEG
          CALL FMDM(X,MA)
          GO TO 150
      ENDIF

!             Normalize MA.

      KB = N1 - J + 2
      JM2 = J - 2
      DO 130 J = 2, KB
         MA(J) = MA(J+JM2)
  130 CONTINUE
      KB1 = KB + 1
      IF (KB1.LE.N1) THEN
          DO 140 J = KB1, N1
             MA(J) = 0
  140     CONTINUE
      ENDIF

      IF (INTEG.LT.0) MA(2) = -MA(2)

  150 RETURN
END SUBROUTINE

SUBROUTINE FMINP(LINE,MA,NPT,LB)

!  Convert an A1 character string to floating point multiple precision
!  format.

!  LINE is an A1 character array of length LB to be converted
!       to FM format and returned in MA.
!  NPT is a pointer telling the routine where in the array to begin
!      the conversion.  This allows more than one number to be stored
!      in an array and converted in place.
!  LB is a pointer to the last character of the field for that number.

!  The input number may be in integer or any real format.
!  The convention is made that if no digits appear before 'E' then 1.0
!  is assumed.  For example 'E6' is converted as '1.0E+6'.

!  So that FMINP will convert any output from FMOUT, LINE is tested
!  to see if the input is one of the special symbols +OVERFLOW,
!  -OVERFLOW, +UNDERFLOW, -UNDERFLOW, or UNKNOWN.
!  For user input the abbreviations OVFL, UNFL, UNKN may be used.

      CHARACTER LINE(LB),KBLANK,KMINUS,KOVFL(4),KUNFL(4),KUNKN(4)
      DIMENSION JTRANS(8,4)
      DIMENSION MA(LUNPCK)

!  Simulate a finite-state automaton to scan the input line
!  and build the number.  States of the machine:

!  1.  Initial entry to the subroutine
!  2.  Sign of the number
!  3.  Scanning digits before a decimal point
!  4.  Decimal point
!  5.  Scanning digits after a decimal point
!  6.  E, D, Q, or M - precision indicator before the exponent
!  7.  Sign of the exponent
!  8.  Scanning exponent
!  9.  Syntax error

!  Character types recognized by the machine:

!  1.  Sign (+,-)
!  2.  Numeral (0,1,...,9)
!  3.  Decimal point (.)
!  4.  Precision indicator (E,D,Q,M)
!  5.  Illegal character for number

!  All blanks are ignored.  The analysis of the number proceeds as
!  follows:  If the simulated machine is in state JSTATE and a character
!  of type JTYPE is encountered the new state of the machine is given by
!  JTRANS(JSTATE,JTYPE).

!  In this DATA statement note the array is loaded by columns.

!          State   1  2  3  4  5  6  7  8
!                                               Type
      DATA JTRANS/ 2, 9, 9, 9, 9, 7, 9, 9,  &
                   3, 3, 3, 5, 5, 8, 8, 8,  &  ! N
                   4, 4, 4, 9, 9, 9, 9, 9,  &  ! D
                   6, 6, 6, 6, 6, 9, 9, 9 /    ! P

      DATA KBLANK/' '/,KMINUS/'-'/,KOVFL/'O','V','F','L'/,  &
           KUNFL/'U','N','F','L'/, KUNKN/'U','N','K','N'/

      NCALL = NCALL + 1
      KSTACK(NCALL) = 17
      IF (ABS(NTRACE).GE.2 .AND. NCALL.LE.LVLTRC) THEN
          WRITE (KW,110) (LINE(L),L=NPT,LB)
  110     FORMAT(' INPUT TO FMINP'/(14X,65A1))
      ENDIF
      NDSAVE = NDIG
      KASAVE = KARGSW
      KARGSW = 1
      KRSAVE = KROUND
      KROUND = 1
      KFLAG = 0

!             Since arithmetic tracing is not usually desired during
!             I/O conversion, disable tracing during this routine.

      NTRSAV = NTRACE
      NTRACE = 0

!             Check for special symbols.

      KMN = 1
      KOF = 1
      KUF = 1
      KUK = 1
      DO 120 J = NPT, LB
         IF (LINE(J).EQ.KMINUS) KMN = -1
         IF (LINE(J).EQ.KOVFL(KOF)) THEN
             KOF = KOF + 1
             IF (KOF.EQ.5) THEN
                 CALL FMIM(0,MA)
                 MA(1) = KEXPOV
                 MA(2) = KMN
                 GO TO 250
             ENDIF
         ENDIF
         IF (LINE(J).EQ.KUNFL(KUF)) THEN
             KUF = KUF + 1
             IF (KUF.EQ.5) THEN
                 CALL FMIM(0,MA)
                 MA(1) = KEXPUN
                 MA(2) = KMN
                 GO TO 250
             ENDIF
         ENDIF
         IF (LINE(J).EQ.KUNKN(KUK)) THEN
             KUK = KUK + 1
             IF (KUK.EQ.5) THEN
                 CALL FMIM(0,MA)
                 MA(1) = KUNKNO
                 MA(2) = 1
                 GO TO 250
             ENDIF
         ENDIF
  120 CONTINUE

!             Increase the working precision.

      IF (NCALL.EQ.1) THEN
          B = JBASE
          K = 5.0*LOG(10.0)/LOG(B) + 2.0
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MA(1) = KUNKNO
              MA(2) = 1
              DO 130 J = 2, NDSAVE
                 MA(J+1) = 0
  130         CONTINUE
              GO TO 250
          ENDIF
      ENDIF
      NDSAV1 = NDIG
      KSTART = NPT
      KSTOP = LB
      JSTATE = 1
      KSIGN = 1
      N1 = NDIG + 1
      DO 140 J = 1, N1
         M02(J) = 0
         M03(J) = 0
         M04(J) = 0
         M05(J) = 0
  140 CONTINUE

!             If JBASE is a power of ten then call FMINP2 for a much
!             faster input conversion.

      KPOWER = LOG10(REAL(JBASE)) + 0.5
      IF (JBASE.EQ.10**KPOWER) THEN
          CALL FMINP2(MA,LINE,KSTART,KSTOP,JTRANS,KBLANK,KPOWER)
          GO TO 240
      ENDIF

      N2 = 0
      KSIGNX = 1
      KF1 = 0
      KF2 = 0
      KEXP = 0
      KTENF1 = 1
      KTENF2 = 1
      KTENEX = 1
      K10PWR = 0

!             Large is a threshold used in order to do as much of the
!             conversion as possible in one-word integer arithmetic.

      LARGE = (MXBASE*MXBASE - 10)/10

!             KDFLAG will be 1 if any digits are found before 'E'.

      KDFLAG = 0

!             Scan the number.

      DO 210 J = KSTART, KSTOP
         IF (LINE(J).EQ.KBLANK) GO TO 210
         CALL FMINPT(LINE(J),KTYPE,KVAL)
         IF (KTYPE.GE.5) GO TO 260

         JSTATE = JTRANS(JSTATE,KTYPE)

         GO TO (260,150,160,210,170,180,190,200,260),JSTATE

!             State 2.  Sign of the number.

  150    KSIGN = KVAL
         GO TO 210

!             State 3.  Digits before a decimal point.

  160    KDFLAG = 1
         KF1 = 10*KF1 + KVAL
         KTENF1 = 10*KTENF1
         IF (KTENF1.GT.LARGE) THEN
             IF (KTENF1.NE.K10PWR .AND. M03(2).NE.0) THEN
                 CALL FMI2M(KTENF1,MA)
                 K10PWR = KTENF1
             ENDIF
             IF (M03(2).EQ.0) THEN
                 CALL FMI2M(KF1,M03)
             ELSE
                 NDIG = MAX(2,MIN(M03(1)+MA(1),NDSAV1))
                 CALL FMMPY(M03,MA,M03)
                 NDIG = NDSAV1
                 CALL FMI2M(KF1,M02)
                 NDIG = MAX(2,MIN(MAX(M03(1),M02(1))+1,NDSAV1))
                 IF (KF1.NE.0) CALL FMADD(M03,M02,M03)
                 NDIG = NDSAV1
             ENDIF
             KF1 = 0
             KTENF1 = 1
         ENDIF
         GO TO 210

!             State 5.  Digits after a decimal point.

  170    KDFLAG = 1
         N2 = N2 + 1
         KF2 = 10*KF2 + KVAL
         KTENF2 = 10*KTENF2
         IF (KTENF2.GT.LARGE) THEN
             IF (KTENF2.NE.K10PWR .AND. M04(2).NE.0) THEN
                 CALL FMI2M(KTENF2,MA)
                 K10PWR = KTENF2
             ENDIF
             IF (M04(2).EQ.0) THEN
                 CALL FMI2M(KF2,M04)
             ELSE
                 NDIG = MAX(2,MIN(M04(1)+MA(1),NDSAV1))
                 CALL FMMPY(M04,MA,M04)
                 NDIG = NDSAV1
                 CALL FMI2M(KF2,M02)
                 NDIG = MAX(2,MIN(MAX(M04(1),M02(1))+1,NDSAV1))
                 IF (KF2.NE.0) CALL FMADD(M04,M02,M04)
                 NDIG = NDSAV1
             ENDIF
             KF2 = 0
             KTENF2 = 1
         ENDIF
         GO TO 210

!             State 6.  Precision indicator.

  180    IF (KDFLAG.EQ.0) CALL FMI2M(1,M03)
         GO TO 210

!             State 7.  Sign of the exponent.

  190    KSIGNX = KVAL
         GO TO 210

!             State 8.  Digits of the exponent.

  200    KEXP = 10*KEXP + KVAL
         KTENEX = 10*KTENEX
         IF (KTENEX.GT.LARGE) THEN
             IF (KTENEX.NE.K10PWR .AND. M05(2).NE.0) THEN
                 CALL FMI2M(KTENEX,MA)
                 K10PWR = KTENEX
             ENDIF
             IF (M05(2).EQ.0) THEN
                 CALL FMI2M(KEXP,M05)
             ELSE
                 NDIG = MAX(2,MIN(M05(1)+MA(1),NDSAV1))
                 CALL FMMPY(M05,MA,M05)
                 NDIG = NDSAV1
                 CALL FMI2M(KEXP,M02)
                 NDIG = MAX(2,MIN(MAX(M05(1),M02(1))+1,NDSAV1))
                 IF (KEXP.NE.0) CALL FMADD(M05,M02,M05)
                 NDIG = NDSAV1
             ENDIF
             KEXP = 0
             KTENEX = 1
         ENDIF

  210 CONTINUE

!             Form the number and return.
!             MA = KSIGN*(M03 + M04/10.0**N2)*10.0**M05

      IF (KTENF1.GT.1) THEN
          IF (KTENF1.NE.K10PWR .AND. M03(2).NE.0) THEN
              CALL FMI2M(KTENF1,MA)
              K10PWR = KTENF1
          ENDIF
          IF (M03(2).EQ.0) THEN
              CALL FMI2M(KF1,M03)
          ELSE
              NDIG = MAX(2,MIN(M03(1)+MA(1),NDSAV1))
              CALL FMMPY(M03,MA,M03)
              NDIG = NDSAV1
              CALL FMI2M(KF1,M02)
              NDIG = MAX(2,MIN(MAX(M03(1),M02(1))+1,NDSAV1))
              IF (KF1.NE.0) CALL FMADD(M03,M02,M03)
              NDIG = NDSAV1
          ENDIF
      ENDIF
      IF (KTENF2.GT.1) THEN
          IF (KTENF2.NE.K10PWR .AND. M04(2).NE.0) THEN
              CALL FMI2M(KTENF2,MA)
              K10PWR = KTENF2
          ENDIF
          IF (M04(2).EQ.0) THEN
              CALL FMI2M(KF2,M04)
          ELSE
              NDIG = MAX(2,MIN(M04(1)+MA(1),NDSAV1))
              CALL FMMPY(M04,MA,M04)
              NDIG = NDSAV1
              CALL FMI2M(KF2,M02)
              NDIG = MAX(2,MIN(MAX(M04(1),M02(1))+1,NDSAV1))
              IF (KF2.NE.0) CALL FMADD(M04,M02,M04)
              NDIG = NDSAV1
          ENDIF
      ENDIF
      IF (KTENEX.GT.1) THEN
          IF (KTENEX.NE.K10PWR .AND. M05(2).NE.0) THEN
              CALL FMI2M(KTENEX,MA)
              K10PWR = KTENEX
          ENDIF
          IF (M05(2).EQ.0) THEN
              CALL FMI2M(KEXP,M05)
          ELSE
              NDIG = MAX(2,MIN(M05(1)+MA(1),NDSAV1))
              CALL FMMPY(M05,MA,M05)
              NDIG = NDSAV1
              CALL FMI2M(KEXP,M02)
              NDIG = MAX(2,MIN(MAX(M05(1),M02(1))+1,NDSAV1))
              IF (KEXP.NE.0) CALL FMADD(M05,M02,M05)
              NDIG = NDSAV1
          ENDIF
      ENDIF

      IF (KSIGNX.EQ.-1) M05(2) = -M05(2)
      IF (M04(2).NE.0) THEN
          CALL FMI2M(10,M02)
          K = N2
          IF (MOD(K,2).EQ.0) THEN
              CALL FMI2M(1,MA)
          ELSE
              CALL FMEQU(M02,MA,NDIG,NDIG)
          ENDIF

  220     K = K/2
          NDIG = MAX(2,MIN(2*M02(1),NDSAV1))
          CALL FMMPY(M02,M02,M02)
          IF (MOD(K,2).EQ.1) THEN
              NDIG = MAX(2,MIN(M02(1)+MA(1),NDSAV1))
              CALL FMMPY(M02,MA,MA)
          ENDIF
          IF (K.GT.1) GO TO 220
          NDIG = NDSAV1
          CALL FMDIV(M04,MA,M04)
      ENDIF
      IF (M05(2).NE.0) THEN
          CALL FMI2M(10,M02)
          KWSV = KWARN
          KWARN = 0
          CALL FMM2I(M05,KEXP)
          KWARN = KWSV
          IF (KFLAG.NE.0) GO TO 260
          K = ABS(KEXP)
          IF (MOD(K,2).EQ.0) THEN
              CALL FMI2M(1,M05)
          ELSE
              CALL FMEQU(M02,M05,NDIG,NDIG)
          ENDIF

  230     K = K/2
          NDIG = MAX(2,MIN(2*M02(1),NDSAV1))
          CALL FMMPY(M02,M02,M02)
          IF (MOD(K,2).EQ.1) THEN
              NDIG = MAX(2,MIN(M02(1)+M05(1),NDSAV1))
              CALL FMMPY(M02,M05,M05)
          ENDIF
          IF (K.GT.1) GO TO 230
          NDIG = NDSAV1
          IF (KEXP.LT.0) THEN
              CALL FMI2M(1,M02)
              CALL FMDIV(M02,M05,M05)
          ENDIF
      ENDIF
      CALL FMADD(M03,M04,MA)
      IF (M05(2).NE.0) CALL FMMPY(MA,M05,MA)
      IF (KSIGN.EQ.-1) MA(2) = -MA(2)
  240 CALL FMEQU(MA,MA,NDIG,NDSAVE)
      IF (MA(1).EQ.KUNKNO) GO TO 260

  250 NDIG = NDSAVE
      KARGSW = KASAVE
      NTRACE = NTRSAV
      KROUND = KRSAVE
      IF (KFLAG.EQ.1) KFLAG = 0
      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
      RETURN

!             Error in converting the number.

  260 CALL FMIM(0,MA)
      MA(1) = KUNKNO
      MA(2) = 1
      KFLAG = -7
      CALL FMWARN
      GO TO 250
END SUBROUTINE

SUBROUTINE FMINP2(MA,LINE,KSTART,KSTOP,JTRANS,KBLANK,KPOWER)

!  Internal routine for input conversion for a power of ten JBASE.

      CHARACTER LINE(KSTOP),KBLANK
      DIMENSION JTRANS(8,4)
      DIMENSION MA(LUNPCK)

      JSTATE = 1
      KDFLAG = 0
      KSIGN = 1
      KSIGNX = 1
      KF1 = 0
      KNZDIG = 0
      KF1DIG = 0
      KF2 = 0
      KF2DIG = 0
      KF2PT = 2
      KEXP = 0
      LARGE = MXBASE*MXBASE/10

!             Scan the number.

      DO 170 J = KSTART, KSTOP
         IF (LINE(J).EQ.KBLANK) GO TO 170
         CALL FMINPT(LINE(J),KTYPE,KVAL)
         IF (KTYPE.GE.5) GO TO 180

         JSTATE = JTRANS(JSTATE,KTYPE)

         GO TO (180,110,120,170,130,140,150,160,180),JSTATE

!             State 2.  Sign of the number.

  110    KSIGN = KVAL
         GO TO 170

!             State 3.  Digits before a decimal point.

  120    KDFLAG = 1
         KF1 = 10*KF1 + KVAL
         IF (KVAL.GT.0 .OR. KNZDIG.NE.0) THEN
             KNZDIG = 1
             KF1DIG = KF1DIG + 1
         ENDIF
         IF (KF1DIG.EQ.KPOWER) THEN
             M03(1) = M03(1) + 1
             IF (M03(1).LT.NDIG) M03(M03(1)+1) = KF1
             KF1 = 0
             KF1DIG = 0
         ENDIF
         GO TO 170

!             State 5.  Digits after a decimal point.

  130    KDFLAG = 1
         IF (KF2PT.GT.NDIG+1) GO TO 170
         KF2 = 10*KF2 + KVAL
         KF2DIG = KF2DIG + 1
         IF (KF2DIG.EQ.KPOWER) THEN
             M04(KF2PT) = KF2
             IF (KF2.EQ.0 .AND. KF2PT.EQ.2) THEN
                 M04(1) = M04(1) - 1
             ELSE
                 KF2PT = KF2PT + 1
             ENDIF
             KF2 = 0
             KF2DIG = 0
         ENDIF
         GO TO 170

!             State 6.  Precision indicator.

  140    IF (KDFLAG.EQ.0) CALL FMI2M(1,M03)
         GO TO 170

!             State 7.  Sign of the exponent.

  150    KSIGNX = KVAL
         GO TO 170

!             State 8.  Digits of the exponent.

  160    IF (KEXP.GE.LARGE) THEN
             IF (M03(2).EQ.0 .AND. M04(2).EQ.0) THEN
                 CALL FMI2M(0,MA)
                 RETURN
             ENDIF
             CALL FMIM(0,MA)
             IF (KSIGNX.EQ.1) THEN
                 MA(1) = KEXPOV
                 KFLAG = -4
             ELSE
                 MA(1) = KEXPUN
                 KFLAG = -4
             ENDIF
             MA(2) = KSIGN
             CALL FMWARN
             RETURN
         ENDIF

         KEXP = 10*KEXP + KVAL

  170 CONTINUE

!             Form the number and return.
!             MA = KSIGN*(M03 + M04)*10.0**(KSIGNX*KEXP)

      IF (KF1DIG.NE.0) THEN
          M03(1) = M03(1) + 1
          KSHIFT = 10**(KPOWER-KF1DIG)
          IF (M03(1).LT.NDIG) M03(M03(1)+1) = KF1*KSHIFT
          IF (KSHIFT.GT.1) CALL FMDIVI(M03,KSHIFT,M03)
      ENDIF

      IF (KF2DIG.NE.0) THEN
          KSHIFT = 10**(KPOWER-KF2DIG)
          M04(KF2PT) = KF2*KSHIFT
      ENDIF
      IF (M04(2).EQ.0) M04(1) = 0

      IF (KEXP.NE.0) THEN
          IF (KSIGNX.EQ.1) THEN
              M05(1) = KEXP/KPOWER + 1
              M05(2) = 10**(MOD(KEXP,KPOWER))
          ELSE
              M05(1) = -((KEXP-1)/KPOWER)
              KSHIFT = 10**(MOD(KEXP,KPOWER))
              IF (KSHIFT.GT.1) THEN
                  M05(2) = JBASE/KSHIFT
              ELSE
                  M05(2) = 1
              ENDIF
          ENDIF
      ENDIF

      CALL FMADD(M03,M04,MA)
      IF (KEXP.GT.0) CALL FMMPY(MA,M05,MA)
      MA(2) = KSIGN*MA(2)

      RETURN

!             Error in converting the number.

  180 CALL FMIM(0,MA)
      MA(1) = KUNKNO
      MA(2) = 1
END SUBROUTINE

SUBROUTINE FMINPT(KCHAR,KTYPE,KVAL)

!  Character KCHAR is looked up in the table of valid characters.
!  KTYPE is returned as its character type.
!  KVAL is returned as the associated value.

!  Types
!    1.  Sign (+,-)
!    2.  Numeral (0,1,...,9)
!    3.  Decimal point (.)
!    4.  Precision indicator (E,D,Q,M)
!    5.  Illegal character for number

      CHARACTER KCHAR
      CHARACTER *21 LCHARS
      DIMENSION LTYPES(21),LVALS(21)
      DATA LCHARS /'+-0123456789.EDQMedqm'/
      DATA LTYPES/1,1,10*2,3,8*4/
      DATA LVALS/1,-1,0,1,2,3,4,5,6,7,8,9,9*0/

      J = INDEX(LCHARS,KCHAR)
      IF (J.GT.0) THEN
          KTYPE = LTYPES(J)
          KVAL = LVALS(J)
      ELSE
          KTYPE = 5
          KVAL = 0
      ENDIF
END SUBROUTINE

SUBROUTINE FMINT(MA,MB)

!  MB = INT(MA)

!  The integer part of MA is computed and returned in MB as a multiple
!  precision floating point number.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 18
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)
      IF (KARGSW.NE.1) THEN
          CALL FMARGS(18,1,MA,MB,KRESLT)
          IF (KRESLT.NE.0) THEN
              CALL FMRSLT(MA,MA,MB,KRESLT)
              IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
              NCALL = NCALL - 1
              RETURN
          ENDIF
      ENDIF

      N1 = NDIG + 1

!             If MA is less than one in magnitude, return zero.

      IF (MA(1).LE.0) THEN
          DO 110 J = 1, N1
             MB(J) = 0
  110     CONTINUE
          GO TO 150
      ENDIF

!             If the radix point is off the right end of MA then MA is
!             already an integer.  Return MA.

      IF (MA(1).GE.NDIG) THEN
          DO 120 J = 1, N1
             MB(J) = MA(J)
  120     CONTINUE
          GO TO 150
      ENDIF

!             Here MA has both integer and fraction parts.  Replace
!             the digits right of the radix point by zeros.

      KA = MA(1) + 2
      KB = KA - 1
      DO 130 J = 1, KB
         MB(J) = MA(J)
  130 CONTINUE

      DO 140 J = KA, N1
         MB(J) = 0
  140 CONTINUE

  150 IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMIPWR(MA,INT_T,MB)

!  MB = MA ** INT_T

!  Raise an FM number to an integer power.
!  The binary multiplication method used requires an average of
!  1.5 * LOG2(INT_T) multiplications.  MA may be negative.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 19
      IF (NTRACE.NE.0) THEN
          CALL FMNTR(2,MA,MA,1)
          CALL FMNTRI(2,INT_T,0)
      ENDIF

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO .OR. (INT_T.LE.0 .AND. MA(2).EQ.0)) THEN
          MA2 = MA(2)
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          IF (INT_T.LE.0 .AND. MA2.EQ.0) CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (INT_T.EQ.0) THEN
          CALL FMIM(1,MB)
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (ABS(INT_T).EQ.1) THEN
          KWRNSV = KWARN
          KWARN = 0
          IF (INT_T.EQ.1) THEN
              CALL FMEQU(MA,MB,NDIG,NDIG)
          ELSE
              CALL FMIM(1,M01)
              CALL FMDIV(M01,MA,MB)
          ENDIF
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          KWARN = KWRNSV
          RETURN
      ENDIF

      IF (MA(2).EQ.0) THEN
          CALL FMIM(0,MB)
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPOV) THEN
          JSIGN = 1
          IF (MA(2).LT.0) JSIGN = -1
          CALL FMIM(0,MB)
          IF (INT_T.GT.0) THEN
              MB(1) = KEXPOV
              MB(2) = JSIGN**MOD(INT_T,2)
              KFLAG = -5
          ELSE
              MB(1) = KEXPUN
              MB(2) = JSIGN**MOD(INT_T,2)
              KFLAG = -6
          ENDIF
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPUN) THEN
          JSIGN = 1
          IF (MA(2).LT.0) JSIGN = -1
          CALL FMIM(0,MB)
          IF (INT_T.GT.0) THEN
              MB(1) = KEXPUN
              MB(2) = JSIGN**MOD(INT_T,2)
              KFLAG = -6
          ELSE
              MB(1) = KEXPOV
              MB(2) = JSIGN**MOD(INT_T,2)
              KFLAG = -5
          ENDIF
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Increase the working precision.

      NDSAVE = NDIG
      IF (NCALL.EQ.1) THEN
          B = JBASE
          XINT = ABS(INT_T)
          K = (5.0*LOG(10.0) + LOG(XINT))/LOG(B) + 2.0
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MB(1) = KUNKNO
              MB(2) = 1
              DO 110 J = 2, NDSAVE
                 MB(J+1) = 0
  110         CONTINUE
              NDIG = NDSAVE
              IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
              NCALL = NCALL - 1
              RETURN
          ENDIF
      ENDIF

!             Initialize.

      KWSV = KWARN
      KWARN = 0
      K = ABS(INT_T)
      CALL FMEQU(MA,M01,NDSAVE,NDIG)
      IF (MOD(K,2).EQ.0) THEN
          CALL FMI2M(1,MB)
      ELSE
          CALL FMEQU(M01,MB,NDIG,NDIG)
      ENDIF

!             This is the multiplication loop.

  120 K = K/2
      CALL FMMPY(M01,M01,M01)
      IF (MOD(K,2).EQ.1) CALL FMMPY(M01,MB,MB)
      IF (K.GT.1) GO TO 120

!             Invert if the exponent is negative.

      IF (INT_T.LT.0) THEN
          CALL FMI2M(1,M01)
          CALL FMDIV(M01,MB,MB)
      ENDIF
      KWARN = KWSV

!             Round the result and return.

      CALL FMEQU(MB,MB,NDIG,NDSAVE)
      NDIG = NDSAVE
      IF (KFLAG.LT.0) CALL FMWARN
      IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMLG10(MA,MB)

!  MB = LOG10(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(20,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,MB,NDSAVE,NDIG)
      CALL FMLN(MB,MB)
      CALL FMLNI(10,M03)
      CALL FMDIV(MB,M03,MB)

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMLN(MA,MB)

!  MB = LOG(MA)     (Natural logarithm)

      DOUBLE PRECISION Y
      DIMENSION MA(LUNPCK),MB(LUNPCK),NSTACK(19)

      CALL FMENTR(21,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

!             If MA is close to 1, use the Taylor series:
!                   LN(1+X) = X - X**2/2 + X**3/3 - ...
!             This is faster for small X and avoids cancellation error.

!             This method is faster for moderate sized NDIG, but is
!             asymptotically slower by a factor of NDIG**(2/3) than
!             using Newton and FMEXP.  For JBASE=10,000 the Taylor
!             series is faster for NDIG less than about 150 (and is
!             used only when MA is between .9999 and 1.0001).

      IF ((MA(1).EQ.1 .AND. MA(2).EQ.1 .AND. MA(3).EQ.0) .OR.  &
          (MA(1).EQ.0 .AND. MA(2).EQ.JBASE-1)) THEN
           CALL FMEQU(MA,M03,NDSAVE,NDIG)
           CALL FMI2M(-1,M01)
           CALL FMADD(M03,M01,M03)

!             The sum will be done as two concurrent series.

           NDSAV1 = NDIG
           CALL FMEQU(M03,M04,NDIG,NDIG)
           CALL FMDIVI(M03,2,M05)
           CALL FMMPY(M03,M03,MB)
           CALL FMEQU(M03,M02,NDIG,NDIG)
           KBOT = 2

  110      KBOT = KBOT + 1
           CALL FMMPY(M02,MB,M02)
           CALL FMDIVI(M02,KBOT,M01)
           NDIG = NDSAV1
           CALL FMADD(M04,M01,M04)
           NDIG = MAX(2,NDSAV1 - (M04(1)-M01(1)))
           KBOT = KBOT + 1
           CALL FMDIVI(M02,KBOT,M01)
           NDIG = NDSAV1
           CALL FMADD(M05,M01,M05)
           NDIG = MAX(2,NDSAV1 - (M04(1)-M01(1)))
           IF (KFLAG.NE.1) GO TO 110

           NDIG = NDSAV1
           CALL FMMPY(M05,M03,M05)
           CALL FMSUB(M04,M05,MB)
           GO TO 160
      ENDIF

!             Check to see if the argument is a small integer.
!             If so use FMLNI.

      KM1 = 0
      MA1 = MA(1)
      CALL FMEQU(MA,M05,NDSAVE,NDIG)
      KWSV = KWARN
      KWARN = 0
      CALL FMM2I(M05,INT_T)
      KWARN = KWSV
      IF (KFLAG.EQ.0 .AND. INT_T.LT.MXBASE) THEN
          CALL FMLNI(INT_T,MB)
          GO TO 160
      ENDIF

!             See if the argument can be scaled to a small integer.

      N3 = NDIG + 3
      N1 = NDIG + 1
      DO 120 J = 2, N1
         IF (M05(N3-J).NE.0) THEN
             LAST = N3 - J - 1
             GO TO 130
         ENDIF
  120 CONTINUE

  130 KSCALE = MA1 - LAST
      M05(1) = LAST
      KWSV = KWARN
      KWARN = 0
      CALL FMM2I(M05,INT_T)
      KWARN = KWSV
      IF (KFLAG.EQ.0 .AND. INT_T.LT.MXBASE) THEN
          CALL FMLNI(INT_T,M04)
          IF (INT_T.EQ.1) KM1 = 1
          K2EXP = 0
          GO TO 150
      ENDIF

!             For the non-integer case scale the argument to lie between
!             e/2 and e to speed up the calls to FMEXP.

      M05(1) = 1
      KSCALE = MA1 - 1
      CALL FMM2DP(M05,Y)
      K2EXP = LOG(2.0*Y/2.71828)/0.693147
      IF (Y.LT.1.359141) THEN
          K2EXP = -1
          CALL FMMPYI(M05,2,M05)
          Y = 2.0*Y
      ELSE
          K2 = 2**K2EXP
          CALL FMDIVI(M05,K2,M05)
          Y = Y/K2
      ENDIF

!             Generate the initial approximation.

      Y = LOG(Y)
      CALL FMDP2M(Y,M04)
      CALL FMI2M(1,MB)
      CALL FMDIG(NSTACK,KST)

!             Newton iteration.

      DO 140 J = 1, KST
         NDIG = NSTACK(J)
         CALL FMEXP(M04,MB)
         CALL FMSUB(M05,MB,M02)
         CALL FMDIV(M02,MB,MB)
         CALL FMADD(M04,MB,M04)
  140 CONTINUE

!             Compute LN(JBASE**KSCALE).

  150 IF ((NJBLB.NE.JBASE .OR. NDIGLB.LT.NDIG) .AND. KSCALE.NE.0) THEN
          NDSV = NDIG
          NDIG = MIN(NDIG+2,MXNDG2)
          CALL FMLNI(JBASE,MLBSAV)
          NJBLB = JBASE
          NDIGLB = NDIG
          NDIG = NDSV
      ENDIF

      IF (KSCALE.NE.0 .AND. KM1.EQ.0) THEN
          CALL FMMPYI(MLBSAV,KSCALE,MB)
          CALL FMADD(M04,MB,MB)
      ELSE IF (KSCALE.NE.0 .AND. KM1.EQ.1) THEN
          CALL FMMPYI(MLBSAV,KSCALE,MB)
      ELSE IF (KSCALE.EQ.0 .AND. KM1.EQ.0) THEN
          CALL FMEQU(M04,MB,NDIG,NDIG)
      ELSE IF (KSCALE.EQ.0 .AND. KM1.EQ.1) THEN
          CALL FMI2M(0,MB)
      ENDIF

      IF (K2EXP.NE.0) THEN
          CALL FMLNI(2,M04)
          CALL FMMPYI(M04,K2EXP,M04)
          CALL FMADD(MB,M04,MB)
      ENDIF

!             Round the result and return.

  160 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMLNI(INT_T,MA)

!  MA = LOG(INT_T)

!  Compute the natural logarithm of an integer INT_T.

!  If INT_T has only powers of 2, 3, 5, and 7 in its factorization then
!  FMLNI is faster than FMLN.  Otherwise, if INT_T.GE.MXBASE (i.e., INT_T
!  does not fit in 1/2 word) then FMLN is usually faster.

!  Use FMLN instead of FMLNI if 10*INT_T would cause integer overflow.

      DIMENSION MA(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 22
      IF (NTRACE.NE.0) CALL FMNTRI(2,INT_T,1)

!             Check for special cases.

      IF (INT_T.LE.0) THEN
          CALL FMIM(0,MA)
          MA(1) = KUNKNO
          MA(2) = 1
          KFLAG = -4
          CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (INT_T.EQ.1) THEN
          CALL FMI2M(0,MA)
          IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Increase the working precision.

      NDSAVE = NDIG
      IF (NCALL.EQ.1) THEN
          B = JBASE
          K = 5.0*LOG(10.0)/LOG(B) + 2.0
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MA(1) = KUNKNO
              MA(2) = 1
              DO 110 J = 2, NDSAVE
                 MA(J+1) = 0
  110         CONTINUE
              NDIG = NDSAVE
              IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
              NCALL = NCALL - 1
              RETURN
          ENDIF
      ENDIF
      KASAVE = KARGSW
      KARGSW = 1

!             Find integers K2, K3, K5, and K7 such that
!                NT = 2**K2 * 3**K3 * 5**K5 * 7**K7
!             is a good approximation of INT_T.
!             KDELTA = ABS(INT_T - NT).

      INT2 = INT_T
      IF (INT_T.GT.MAXINT/100) INT2 = INT_T/100
      KDELTA = INT2
      NT = 0
      K2 = 0
      K3 = 0
      K5 = 0
      K7 = 0

!             Start the search loop.

      XINT = INT2
      LAST = LOG(XINT)/LOG(2.0) + 2.0

      JTEMP7 = 1
      DO 180 J7 = 1, LAST
         IF (JTEMP7.GT.INT2 .AND.  &
             ABS(JTEMP7-INT2).GT.KDELTA) GO TO 190

         JTEMP5 = JTEMP7
         DO 160 J5 = 1, LAST
            IF (JTEMP5.GT.INT2 .AND.  &
                ABS(JTEMP5-INT2).GT.KDELTA) GO TO 170

            JTEMP3 = JTEMP5
            DO 140 J3 = 1, LAST
               IF (JTEMP3.GT.INT2 .AND.  &
                   ABS(JTEMP3-INT2).GT.KDELTA) GO TO 150

               JTEMP2 = JTEMP3
               DO 120 J2 = 1, LAST
                  IF (ABS(JTEMP2-INT2).LE.KDELTA) THEN
                      IF (ABS(JTEMP2-INT2).EQ.KDELTA .AND.  &
                          JTEMP2.LT.INT2)GO TO 130
                      KDELTA = ABS(JTEMP2-INT2)
                      NT = JTEMP2
                      K2 = J2 - 1
                      K3 = J3 - 1
                      K5 = J5 - 1
                      K7 = J7 - 1
                      IF (KDELTA.EQ.0) GO TO 190
                  ENDIF
                  IF (JTEMP2.GT.INT2) GO TO 130

                  JTEMP2 = 2*JTEMP2
  120          CONTINUE

  130          JTEMP3 = 3*JTEMP3
  140       CONTINUE

  150       JTEMP5 = 5*JTEMP5
  160    CONTINUE

  170    JTEMP7 = 7*JTEMP7
  180 CONTINUE

!             If INT_T was too close to the integer overflow limit,
!             restore NT to an approximation of INT_T.

  190 IF (INT2.NE.INT_T) THEN
          IF (NT.LE.INT2) THEN
              NT = NT*100
              K2 = K2 + 2
              K5 = K5 + 2
          ELSE IF (NT.LE.INT_T/98) THEN
              NT = NT*98
              K2 = K2 + 1
              K7 = K7 + 2
          ELSE
              NT = NT*70
              K2 = K2 + 1
              K5 = K5 + 1
              K7 = K7 + 1
          ENDIF
      ENDIF

!             End of the search.  Now compute LN(NT) as a linear
!             combination of LN(125/126), LN(224/225), LN(2400/2401),
!             and LN(4374/4375).  Each of these arguments is near 1.0
!             and each logarithm is computed using only addition and
!             integer division operations -- both O(NDIG) speed.

      KA = -72*K2 - 114*K3 - 167*K5 - 202*K7
      KB = -27*K2 -  43*K3 -  63*K5 -  76*K7
      KC =  19*K2 +  30*K3 +  44*K5 +  53*K7
      KD = -31*K2 -  49*K3 -  72*K5 -  87*K7

      IF (JBASE.NE.NJBLI .OR. NDIG.GT.NDIGLI) THEN
          NDSV = NDIG
          NDIG = MIN(NDIG+2,MXNDG2)
          NJBLI = JBASE
          NDIGLI = NDIG
          CALL FMLNI2(1,126,MLN1)
          CALL FMLNI2(1,225,MLN2)
          CALL FMLNI2(1,2401,MLN3)
          CALL FMLNI2(1,4375,MLN4)
          NDIG = NDSV
      ENDIF

!             If NT.NE.INT_T then the final step is to compute LN(INT_T/NT)
!             and then use LN(INT_T) = LN(INT_T/NT) + LN(NT).

      IF (NT.NE.INT_T) THEN
          ND = NT - INT_T
          CALL FMLNI2(ND,NT,MA)
      ENDIF

      CALL FMMPYI(MLN1,KA,M02)
      CALL FMMPYI(MLN2,KB,M01)
      CALL FMADD(M02,M01,M02)
      CALL FMMPYI(MLN3,KC,M01)
      CALL FMADD(M02,M01,M02)
      CALL FMMPYI(MLN4,KD,M01)
      IF (NT.NE.INT_T) CALL FMADD(M02,MA,M02)
      CALL FMADD(M02,M01,MA)

!             Round and move the result to MA.

      CALL FMEQU(MA,MA,NDIG,NDSAVE)
      NDIG = NDSAVE
      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
      KARGSW = KASAVE
END SUBROUTINE

SUBROUTINE FMLNI2(INT1,INT2,MA)

!  MA = LN(1 - INT1/INT2)

!  Taylor series for computing the logarithm of a rational number
!  near 1.

      DIMENSION MA(LUNPCK)

      CALL FMI2M(INT1,M02)
      CALL FMDIVI(M02,INT2,M02)
      CALL FMEQU(M02,MA,NDIG,NDIG)
      NDSAVE = NDIG
      J = 1

  110 J = J + 1
      IF (INT1.NE.1) CALL FMMPYI(M02,INT1,M02)
      CALL FMDIVI(M02,INT2,M02)
      CALL FMDIVI(M02,J,M01)
      NDIG = NDSAVE
      CALL FMADD(MA,M01,MA)
      NDIG = NDSAVE - (MA(1)-M01(1))
      IF (NDIG.LT.2) NDIG = 2
      IF (KFLAG.NE.1) GO TO 110

      MA(2) = -MA(2)
      NDIG = NDSAVE
END SUBROUTINE

SUBROUTINE FMM2DP(MA,X)

!  X = MA

!  Convert an FM number to double precision.

!  If KFLAG = -4 is returned for a value of MA which is in the range
!  of the machine's double precision number system, change the
!  definition of DPMAX in routine FMSET to reflect the current machine's
!  range.

      DOUBLE PRECISION X
      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 23
      CALL FMARGS(23,1,MA,MA,KRESLT)
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)
      IF (KRESLT.NE.0) THEN

!             Here no valid result can be returned.  Set X to some
!             value which the user is likely to recognize as wrong.

          X = RUNKNO
          KFLAG = -4
          IF (MA(1).NE.KUNKNO) CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTRR(1,X,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      CALL FMMD(MA,X)

      IF (NTRACE.NE.0) CALL FMNTRR(1,X,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMM2I(MA,INT_T)

!  INT_T = MA

!  Convert an FM number to integer.

!  KFLAG =  0 is returned if the conversion is exact.
!        = -4 is returned if MA is larger than MAXINT in magnitude.
!             INT_T = IUNKNO is returned as an indication that INT_T
!             could not be computed without integer overflow.
!        =  2 is returned if MA is smaller than MXBASE**2 in magnitude
!             but MA is not an integer.  The next integer toward zero
!             is returned in INT_T.
!  It is sometimes convenient to call FMM2I in order to see if an FM
!  number can be represented as a one-word integer by checking KFLAG
!  upon return.  To avoid an unwanted error message being printed in
!  the KFLAG=-4 case, set KWARN=0 before the call to FMM2I and reset
!  it after the call.

!  This routine performs the trace printing for the conversion.
!  FMMI is used to do the arithmetic.

      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 24
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)

      CALL FMMI(MA,INT_T)

      IF (NTRACE.NE.0) CALL FMNTRI(1,INT_T,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMM2SP(MA,X)

!  X = MA

!  Convert an FM number to single precision.

!  MA is converted and the result is returned in X.

!  If KFLAG = -4 is returned for a value of MA which is in the range
!  of the machine's single precision number system, change the
!  definition of SPMAX in routine FMSET to reflect the current machine's
!  range.

      DOUBLE PRECISION Y
      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 25
      CALL FMARGS(25,1,MA,MA,KRESLT)
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)
      IF (KRESLT.NE.0) THEN

!             Here no valid result can be returned.  Set X to some
!             value which the user is likely to recognize as wrong.

          X = RUNKNO
          KFLAG = -4
          IF (MA(1).NE.KUNKNO) CALL FMWARN
          Y = X
          IF (NTRACE.NE.0) CALL FMNTRR(1,Y,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      CALL FMMD(MA,Y)
      X = Y

      Y = X
      IF (NTRACE.NE.0) CALL FMNTRR(1,Y,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMMAX(MA,MB,MC)

!  MC = MAX(MA,MB)

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 26
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      KWSV = KWARN
      KWARN = 0
      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
      ELSE IF (FMCOMP(MA,'LT',MB)) THEN
          CALL FMEQU(MB,MC,NDIG,NDIG)
      ELSE
          CALL FMEQU(MA,MC,NDIG,NDIG)
      ENDIF

      KWARN = KWSV
      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMMD(MA,X)

!  X = MA

!  Internal routine for conversion to double precision.

      DOUBLE PRECISION X,Y,YT,XBASE,RZERO,ONE,PMAX
      DIMENSION MA(LUNPCK)

!             Check to see if MA is in range for single or double
!             precision.

      PMAX = DPMAX
      IF (NCALL.GT.0) THEN
          IF (KSTACK(NCALL).EQ.25) PMAX = SPMAX
      ENDIF
      DLOGJB = LOG(DBLE(JBASE))
      DLOGDP = LOG(PMAX)
      MA1 = MA(1)
      NCASE = 0
      IF (DBLE(MA(1)-1)*DLOGJB.GT.DLOGDP) THEN
          KFLAG = -4
          X = RUNKNO
          CALL FMWARN
          RETURN
      ELSE IF (DBLE(MA(1)+1)*DLOGJB.GT.DLOGDP) THEN
          MA(1) = MA(1) - 2
          NCASE = 1
      ELSE IF (DBLE(MA(1)+1)*DLOGJB.LT.-DLOGDP) THEN
          KFLAG = -10
          X = 0.0D0
          CALL FMWARN
          RETURN
      ELSE IF (DBLE(MA(1)-1)*DLOGJB.LT.-DLOGDP) THEN
          MA(1) = MA(1) + 2
          NCASE = 2
      ENDIF

!             Try FMMI first so that small integers will be
!             converted exactly.

      KWSV = KWARN
      KWARN = 0
      CALL FMMI(MA,J)
      KWARN = KWSV
      IF (KFLAG.EQ.0) THEN
          X = J
          RETURN
      ENDIF
      KFLAG = 0

      MA2 = MA(2)
      MA(2) = ABS(MA2)
      RZERO = 0.0
      ONE = 1.0
      N1 = NDIG + 1
      XBASE = JBASE
      X = RZERO
      Y = ONE
      DO 110 J = 2, N1
         Y = Y/XBASE
         YT = MA(J)
         X = X + Y*YT
         YT = ONE + Y*XBASE
         IF (YT.LE.ONE) GO TO 120
  110 CONTINUE

  120 X = X*XBASE**MA(1)
      IF (MA2.LT.0) X = -X
      MA(2) = MA2

!             Check the result if it is near overflow or underflow.

      IF (NCASE.EQ.1) THEN
          IF (X.LE.PMAX/(XBASE*XBASE)) THEN
              X = X*XBASE*XBASE
          ELSE
              KFLAG = -4
              X = RUNKNO
              CALL FMWARN
          ENDIF
      ELSE IF (NCASE.EQ.2) THEN
          IF (X.GE.(1.0D0/PMAX)*XBASE*XBASE) THEN
              X = X/(XBASE*XBASE)
          ELSE
              KFLAG = -10
              X = 0.0D0
              CALL FMWARN
          ENDIF
      ENDIF
      MA(1) = MA1
END SUBROUTINE

SUBROUTINE FMMI(MA,INT_T)

!  INT_T = MA.  Internal FM to integer conversion routine.

      DIMENSION MA(LUNPCK)

      KFLAG = 0
      N1 = NDIG + 1
      LARGE = MAXINT/JBASE
      INT_T = 0
      IF (MA(1).LE.0) THEN
          IF (MA(2).NE.0) KFLAG = 2
          RETURN
      ENDIF

      KB = MA(1) + 1
      INT_T = ABS(MA(2))
      IF (KB.GE.3) THEN
          DO 110 J = 3, KB
             IF (INT_T.GT.LARGE) THEN
                 KFLAG = -4
                 IF (MA(1).NE.KUNKNO) CALL FMWARN
                 INT_T = IUNKNO
                 RETURN
             ENDIF
             IF (J.LE.N1) THEN
                 INT_T = INT_T*JBASE
                 IF (INT_T.GT.MAXINT-MA(J)) THEN
                     KFLAG = -4
                     IF (MA(1).NE.KUNKNO) CALL FMWARN
                     INT_T = IUNKNO
                     RETURN
                 ELSE
                     INT_T = INT_T + MA(J)
                 ENDIF
             ELSE
                 INT_T = INT_T*JBASE
             ENDIF
  110     CONTINUE
      ENDIF

      IF (MA(2).LT.0) INT_T = -INT_T

!             Check to see if MA is an integer.

      KA = KB + 1
      IF (KA.LE.N1) THEN
          DO 120 J = KA, N1
             IF (MA(J).NE.0) THEN
                 KFLAG = 2
                 RETURN
             ENDIF
  120     CONTINUE
      ENDIF

END SUBROUTINE

SUBROUTINE FMMIN(MA,MB,MC)

!  MC = MIN(MA,MB)

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 27
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      KWSV = KWARN
      KWARN = 0
      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
      ELSE IF (FMCOMP(MA,'GT',MB)) THEN
          CALL FMEQU(MB,MC,NDIG,NDIG)
      ELSE
          CALL FMEQU(MA,MC,NDIG,NDIG)
      ENDIF

      KWARN = KWSV
      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMMOD(MA,MB,MC)

!  MC = MA(MOD MB).

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      CALL FMENTR(28,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN
      KARGSW = 0
      KWSV = KWARN
      KWARN = 0

      IF (MB(1).GT.MA(1) .AND. MB(2).NE.0) THEN
          CALL FMEQU(MA,M01,NDIG,NDSAVE)
      ELSE
          IF (NCALL.LE.1 .AND. MA(2).NE.0) NDIG = NDIG + MA(1) - MB(1)
          IF (NDIG.GT.MXNDG2 .OR. MB(2).EQ.0) THEN
              KFLAG = -9
              IF (MA(1).EQ.KEXPOV .OR. MB(1).EQ.KEXPUN .OR. MB(2).EQ.0)  &
                  KFLAG = -4
              KWARN = KWSV
              KARGSW = KASAVE
              MXEXP = MXSAVE
              CALL FMWARN
              MC(1) = KUNKNO
              MC(2) = 1
              DO 110 J = 2, NDSAVE
                 MC(J+1) = 0
  110         CONTINUE
              NDIG = NDSAVE
              IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
              NCALL = NCALL - 1
              RETURN
          ENDIF

          CALL FMEQU(MA,M02,NDSAVE,NDIG)
          CALL FMEQU(MB,M03,NDSAVE,NDIG)
          M02(2) = ABS(M02(2))
          M03(2) = ABS(M03(2))
          CALL FMDIV(M02,M03,M01)
          CALL FMINT(M01,M01)
          CALL FMMPY(M01,M03,M01)
          CALL FMSUB(M02,M01,M01)

!             Due to rounding, M01 may not be between 0 and MB here.

          NTRSAV = NTRACE
          NTRACE = 0
          IF (FMCOMP(M01,'GE',M03)) THEN
              NTRACE = NTRSAV
              CALL FMSUB(M01,M03,M01)
          ENDIF
          NTRACE = NTRSAV
          IF (M01(2).LT.0) CALL FMADD(M01,M03,M01)
          IF (MA(2).LT.0 .AND. M01(1).NE.KUNKNO) M01(2) = -M01(2)
      ENDIF

      IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(M01,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMMOVE(MA)

!  Move a result from the work area (MWA) to MA.

!  If the result has MWA(2)=0, then it is shifted and the exponent
!  adjusted when it is moved to MA.

      DIMENSION MA(LUNPCK)

      IF (MWA(2).NE.0) THEN
          N1 = NDIG + 1
          DO 110 J = 1, N1
             MA(J) = MWA(J)
  110     CONTINUE
      ELSE
          N2 = NDIG + 2
          DO 120 J = 3, N2
             MA(J-1) = MWA(J)
  120     CONTINUE
          IF (MA(2).NE.0) THEN
              MA(1) = MWA(1) - 1
          ELSE
              MA(1) = 0
          ENDIF
      ENDIF

      IF (ABS(MA(1)).GT.MXEXP) CALL FMTRAP(MA)

END SUBROUTINE

SUBROUTINE FMMPY(MA,MB,MC)

!  MC = MA * MB

!  When one of the numbers MA, MB is known to have more zero digits
!  (base JBASE) than the other, it is faster if MB is the one with
!  more zero digits.

!  This routine performs the trace printing for multiplication.
!  FMMPY2 is used to do the arithmetic.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 29
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      CALL FMMPY2(MA,MB,MC)

      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMMPY2(MA,MB,MC)

!  Internal multiplication routine.  MC = MA * MB

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      IF (KARGSW.NE.1) THEN
          CALL FMARGS(29,2,MA,MB,KRESLT)
          IF (KRESLT.NE.0) THEN
              CALL FMRSLT(MA,MB,MC,KRESLT)
              RETURN
          ENDIF
      ELSE
          IF (MA(2).EQ.0 .OR. MB(2).EQ.0) THEN
              CALL FMIM(0,MC)
              RETURN
          ENDIF
      ENDIF
      KFLAG = 0

!             NGUARD is the number of guard digits used.

      IF (NCALL.GT.1) THEN
          NGUARD = 1
      ELSE
          B = JBASE
          NGUARD = 5.0*LOG(10.0)/LOG(B) + 2.0
          IF (NGUARD.GT.NDIG) NGUARD = NDIG
      ENDIF

!             Save the sign of MA and MB and then work only with
!             positive numbers.

      MA2 = MA(2)
      MB2 = MB(2)
      MA(2) = ABS(MA(2))
      MB(2) = ABS(MB(2))

      N1 = NDIG + 1

!             If there is a good chance of finding several zero digits,
!             see which number has more zero digits.

      NZMA = 0
      NZMB = 0
      IF (NDIG/JBASE.GE.6) THEN
          DO 110 J = 2, N1
             IF (MA(J).EQ.0) NZMA = NZMA + 1
             IF (MB(J).EQ.0) NZMB = NZMB + 1
  110     CONTINUE
      ENDIF

!             It is faster if the second argument is the one with
!             more zero digits.

      IF (NZMA.GT.NZMB) THEN
          CALL FMMPY3(MB,MA,NGUARD,KSHIFT)
      ELSE
          CALL FMMPY3(MA,MB,NGUARD,KSHIFT)
      ENDIF

!             The multiplication is complete.  Round the result,
!             move it to MC, and append the correct sign.

      MA(2) = MA2
      MB(2) = MB2
      CALL FMRND(NDIG,NGUARD,KSHIFT)
      CALL FMMOVE(MC)
      IF (MA2*MB2.LT.0) MC(2) = -MC(2)
END SUBROUTINE

SUBROUTINE FMMPY3(MA,MB,NGUARD,KSHIFT)

!  Internal multiplication of MA*MB.  The result is returned in MWA.
!  Both MA and MB are positive.

!  NGUARD is the number of guard digits which will be used.
!  KSHIFT = 1 is returned if a left shift is pending (i.e., MWA(2)=0).
!             The shift will be done in FMMOVE.  KSHIFT = 0 is returned
!             if no shift is pending.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      N1 = NDIG + 1
      MWA(1) = MA(1) + MB(1)
      L = NDIG + 1 + NGUARD
      MWA(L+1) = 0

!             The multiplication loop begins here.

!             NBNORM is the minimum number of digits which can be
!                    multiplied before normalization is required.
!             MAXMWA is an upper bound on the size of values in MWA
!                    divided by (JBASE-1).  It is used to determine
!                    whether to normalize before the next digit is
!                    multiplied.

      NBNORM = MXBASE**2/(JBASE-1)**2
      JBM1 = JBASE - 1
      LMAX = (MXBASE*MXBASE)/JBM1 - JBM1
      IF (NBNORM.GT.1) THEN
          MBJ = MB(2)

!             Count the trailing zeros in MA.

          KNZ = N1
          N3 = NDIG + 3
          DO 110 J = 2, N1
             IF (MA(N3-J).NE.0) THEN
                 KNZ = N3 - J
                 GO TO 120
             ENDIF
  110     CONTINUE

  120     MWA(2) = 0
          N2 = NDIG + 2
          DO 130 K = N2, L
             MWA(K) = 0
  130     CONTINUE
          DO 140 K = 2, N1
             MWA(K+1) = MA(K)*MBJ
  140     CONTINUE
          MAXMWA = MBJ
          DO 170 J = 3, N1
             MBJ = MB(J)
             IF (MBJ.NE.0) THEN
                 MAXMWA = MAXMWA + MBJ
                 JM1 = J - 1
                 KL = MIN(KNZ,L-JM1)

!                       This is the inner loop for multiplication.
!                       KL.GE.2 since NGUARD.GE.1.

                 DO 150 K = 2, KL
                    MWA(JM1+K) = MWA(JM1+K) + MA(K)*MBJ
  150            CONTINUE
             ENDIF

             IF (MAXMWA.GT.LMAX) THEN
                 MAXMWA = 0

!                       Here normalization is only required for the
!                       range of digits currently changing in MWA.

                 JA = JM1 + 2
                 JB = JM1 + KL
                 KB = JB + 1
                 DO 160 K = JA, JB
                    KB = KB - 1
                    MWA(KB-1) = MWA(KB-1) + MWA(KB)/JBASE
                    MWA(KB) = MOD(MWA(KB),JBASE)
  160            CONTINUE
             ENDIF
  170     CONTINUE

!             Perform the final normalization.

          KB = L + 1
          DO 180 K = 3, L
             KB = KB - 1
             MWA(KB-1) = MWA(KB-1) + MWA(KB)/JBASE
             MWA(KB) = MOD(MWA(KB),JBASE)
  180     CONTINUE

      ELSE

!             If normalization must be done for each digit, combine
!             the two loops and normalize as the digits are multiplied.

          DO 190 J = 2, L
             MWA(J) = 0
  190     CONTINUE
          KJ = NDIG + 2
          DO 210 J = 2, N1
             KJ = KJ - 1
             MBKJ = MB(KJ)
             IF (MBKJ.EQ.0) GO TO 210
             KL = L - KJ + 1
             IF (KL.GT.N1) KL = N1
             KI = KL + 2
             KWA = KL+ KJ + 1
             KK = 0
             DO 200 K = 2, KL
                KT = MA(KI-K)*MBKJ + MWA(KWA-K) + KK
                KK = KT/JBASE
                MWA(KWA-K) = KT - JBASE*KK
  200        CONTINUE
             MWA(KWA-KL-1) = KK
  210     CONTINUE

      ENDIF

!             Set KSHIFT = 1 if a shift left is necessary.

      KSHIFT = 0
      IF (MWA(2).EQ.0) KSHIFT = 1

END SUBROUTINE

SUBROUTINE FMMPYI(MA,INT_T,MB)

!  MB = MA * INT_T

!  Multiply FM number MA by one word integer INT_T.

!  This routine is faster than FMMPY when INT_T*JBASE is a
!  one word integer.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 30
      IF (NTRACE.NE.0) THEN
          CALL FMNTR(2,MA,MA,1)
          CALL FMNTRI(2,INT_T,0)
      ENDIF
      KFLAG = 0
      N1 = NDIG + 1

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO) THEN
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (INT_T.EQ.0) THEN
          CALL FMIM(0,MB)
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (ABS(INT_T).EQ.1) THEN
          DO 110 J = 1, N1
             MB(J) = MA(J)
  110     CONTINUE
          IF (MA(1).EQ.KEXPOV) KFLAG = -5
          IF (MA(1).EQ.KEXPUN) KFLAG = -6
          MB(2) = MA(2)*INT_T
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPOV) THEN
          MA2 = MA(2)
          CALL FMIM(0,MB)
          KFLAG = -5
          MB(1) = KEXPOV
          MB(2) = 1
          IF ((MA2.LT.0 .AND. INT_T.GT.0) .OR.  &
              (MA2.GT.0 .AND. INT_T.LT.0))      MB(2) = -1
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPUN) THEN
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          CALL FMWARN
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

!             Work with positive numbers.

      MA2 = MA(2)
      MA(2) = ABS(MA(2))
      INTP = ABS(INT_T)

!             To leave room for the normalization, shift the product
!             to the right KSHIFT places in MWA.

      KSHIFT = (LOG(REAL(MA(2)+1)) + LOG(REAL(INTP))) /  &
                LOG(REAL(JBASE))

!             If INTP is too big use FMMPY.

      IF (KSHIFT.GT.NDIG .OR. INTP.GT.MXBASE**2/JBASE) THEN
          CALL FMIM(INT_T,M01)
          MA(2) = MA2
          CALL FMMPY2(MA,M01,MB)
          IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
          NCALL = NCALL - 1
          RETURN
      ENDIF

      MWA(1) = MA(1) + KSHIFT
      KA = 2 + KSHIFT
      KB = NDIG + 1 + KSHIFT
      KC = NDIG + 5
      IF (KB.LE.KC) THEN
          DO 120 J = KB, KC
             MWA(J) = 0
  120     CONTINUE
      ENDIF

      KCARRY = 0
      KL = KB + KA
      KM = KL - KSHIFT

!             This is the main multiplication loop.

      DO 130 J = KA, KB
         KT = MA(KM-J)*INTP + KCARRY
         MWA(KL-J) = MOD(KT,JBASE)
         KCARRY = KT/JBASE
  130 CONTINUE

!             Resolve the final carry.

      KB = KA - 1
      KL = KB + 2
      IF (2.LE.KB) THEN
          DO 140 J = 2, KB
             MWA(KL-J) = MOD(KCARRY,JBASE)
             KCARRY = KCARRY/JBASE
  140     CONTINUE
      ENDIF

!             Now the first significant digit in the product is in
!             MWA(2) or MWA(3).  Round the result and move it to MB.

      MA(2) = MA2
      IF (MWA(2).EQ.0) THEN
          NGUARD = KSHIFT - 1
          CALL FMRND(NDIG,NGUARD,1)
      ELSE
          CALL FMRND(NDIG,KSHIFT,0)
      ENDIF
      CALL FMMOVE(MB)

!             Put the sign on the result.

      IF ((INT_T.GT.0 .AND. MA2.LT.0) .OR. (INT_T.LT.0 .AND.MA2.GT.0))  &
           MB(2) = -MB(2)

      IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMNAME(NAME)

!  Convert the FM routine number in KSTACK(NCALL) to an A4 character
!  string for output.

      CHARACTER *4 NAME,NROUT
      DIMENSION NROUT(43)

      DATA NROUT/ 'ABS ',  'ACOS',  'ADD ',  'ASIN',  'ATAN',  'ATN2',  &
                  'BIG ',  'COMP',  'COS ',  'COSH',  'DIM ',  'DIV ',  &
                  'DIVI',  'DP2M',  'EXP ',  'I2M ',  'INP ',  'INT ',  &
                  'IPWR',  'LG10',  'LN  ',  'LNI ',  'M2DP',  'M2I ',  &
                  'M2SP',  'MAX ',  'MIN ',  'MOD ',  'MPY ',  'MPYI',  &
                  'NINT',  'OUT ',  'PI  ',  'PWR ',  'SIGN',  'SIN ',  &
                  'SINH',  'SP2M',  'SQRT',  'SUB ',  'TAN ',  'TANH',  &
                  'ULP '  /

      K = KSTACK(NCALL)
      NAME = NROUT(K)
END SUBROUTINE

SUBROUTINE FMNINT(MA,MB)

!  MB = NINT(MA)  --  MB is returned as the nearest integer to MA.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      CALL FMENTR(31,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      KWSV = KWARN
      KWARN = 0
      CALL FMEQU(MA,MB,NDSAVE,NDIG)
      MA2 = MA(2)
      MB(2) = ABS(MB(2))
      CALL FMI2M(1,M01)
      CALL FMDIVI(M01,2,M01)
      CALL FMADD(MB,M01,MB)
      CALL FMINT(MB,MB)
      IF (MA2.LT.0) MB(2) = -MB(2)
      KWARN = KWSV

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMNTR(NTR,MA,MB,NARG)

!  Print FM numbers in base 10 format using FMOUT for conversion.
!  This is used for trace output from the FM routines.

!  NTR =  1 if a result of an FM call is to be printed.
!      =  2 to print input argument(s) to an FM call.

!  MA  -  the FM number to be printed.

!  MB  -  an optional second FM number to be printed.

!  NARG - the number of arguments.  NARG = 1 if only MA is to be
!         printed, and NARG = 2 if both MA and MB are to be printed.


!  NTRACE and LVLTRC (in COMMON /FMUSER/) control trace printout.

!  NTRACE = 0        No printout except warnings and errors.

!  NTRACE = 1        The result of each call to one of the routines
!                    is printed in base 10, using FMOUT.

!  NTRACE = -1       The result of each call to one of the routines
!                    is printed in internal base JBASE format.

!  NTRACE = 2        The input arguments and result of each call to one
!                    of the routines is printed in base 10, using FMOUT.

!  NTRACE = -2       The input arguments and result of each call to one
!                    of the routines is printed in base JBASE format.

!  LVLTRC defines the call level to which the trace is done.  LVLTRC = 1
!         means only FM routines called directly by the user are traced,
!         LVLTRC = 2 also prints traces for FM routines called by other
!         FM routines called directly by the user, etc.

      CHARACTER *4 NAME
      DIMENSION MA(LUNPCK),MB(LUNPCK)

      IF (NTRACE.EQ.0) RETURN
      IF (NCALL.GT.LVLTRC) RETURN
      IF (NTR.EQ.2 .AND. ABS(NTRACE).EQ.1) RETURN

      IF (NTR.EQ.2) THEN
          CALL FMNAME(NAME)
          WRITE (KW,110) NAME
  110     FORMAT(' INPUT TO FM',A4)
      ELSE
          CALL FMNAME(NAME)
          IF (KFLAG.EQ.0) THEN
              WRITE (KW,120) NAME,NCALL,JBASE,NDIG
  120         FORMAT(' FM',A4,15X,'CALL LEVEL =',I2,5X,'JBASE =',  &
                     I10,5X,'NDIG =',I6)
          ELSE
              WRITE (KW,130) NAME,NCALL,JBASE,NDIG,KFLAG
  130         FORMAT(' FM',A4,6X,'CALL LEVEL =',I2,4X,'JBASE =',  &
                     I10,4X,'NDIG =',I6,4X,'KFLAG =',I3)
          ENDIF
      ENDIF

!             Check for base JBASE internal format trace.

      IF (NTRACE.LT.0) THEN
          CALL FMNTRJ(MA,NDIG)
          IF (NARG.EQ.2) CALL FMNTRJ(MB,NDIG)
      ENDIF

!             Check for base 10 trace using FMOUT.

      IF (NTRACE.GT.0) THEN
          CALL FMPRNT(MA)
          IF (NARG.EQ.2) CALL FMPRNT(MB)
      ENDIF

END SUBROUTINE

SUBROUTINE FMNTRI(NTR,N,KNAM)

!  Internal routine for trace output of integer variables.

!  NTR = 1 for output values
!        2 for input values

!  N     Integer to be printed.

!  KNAM  is positive if the routine name is to be printed.

      CHARACTER *4 NAME

      IF (NTRACE.EQ.0) RETURN
      IF (NCALL.GT.LVLTRC) RETURN
      IF (NTR.EQ.2 .AND. ABS(NTRACE).EQ.1) RETURN

      IF (NTR.EQ.2 .AND. KNAM.GT.0) THEN
          CALL FMNAME(NAME)
          WRITE (KW,110) NAME
  110     FORMAT(' INPUT TO FM',A4)
      ENDIF
      IF (NTR.EQ.1 .AND. KNAM.GT.0) THEN
          CALL FMNAME(NAME)
          IF (KFLAG.EQ.0) THEN
              WRITE (KW,120) NAME,NCALL,JBASE,NDIG
  120         FORMAT(' FM',A4,15X,'CALL LEVEL =',I2,5X,'JBASE =',  &
                     I10,5X,'NDIG =',I6)
          ELSE
              WRITE (KW,130) NAME,NCALL,JBASE,NDIG,KFLAG
  130         FORMAT(' FM',A4,6X,'CALL LEVEL =',I2,4X,'JBASE =',  &
                     I10,4X,'NDIG =',I6,4X,'KFLAG =',I3)
          ENDIF
      ENDIF

      WRITE (KW,140) N
  140 FORMAT(1X,I18)

END SUBROUTINE

SUBROUTINE FMNTRJ(MA,ND)

!  Print trace output in internal base JBASE format.  The number to
!  be printed is in MA.

!  ND is the number of base JBASE digits to be printed.

      DIMENSION MA(LUNPCK)

      N1 = ND + 1

      IF (JBASE.LE.10) THEN
          IF (ND.LE.25) THEN
              WRITE (KW,110) (MA(J),J=1,N1)
  110         FORMAT(1X,I19,I4,24I2)
          ELSE
              WRITE (KW,120) (MA(J),J=1,N1)
  120         FORMAT(1X,I19,I4,24I2/(22X,25I2))
          ENDIF
          RETURN
      ENDIF

      IF (JBASE.LE.100) THEN
          IF (ND.LE.15) THEN
              WRITE (KW,130) (MA(J),J=1,N1)
  130         FORMAT(1X,I19,I5,14I3)
          ELSE
              WRITE (KW,140) (MA(J),J=1,N1)
  140         FORMAT(1X,I19,I5,14I3/(22X,15I3))
          ENDIF
          RETURN
      ENDIF

      IF (JBASE.LE.1000) THEN
          IF (ND.LE.10) THEN
              WRITE (KW,150) (MA(J),J=1,N1)
  150         FORMAT(1X,I19,I6,9I4)
          ELSE
              WRITE (KW,160) (MA(J),J=1,N1)
  160         FORMAT(1X,I19,I6,9I4/(22X,10I4))
          ENDIF
          RETURN
      ENDIF

      IF (JBASE.LE.10000) THEN
          IF (ND.LE.10) THEN
              WRITE (KW,170) (MA(J),J=1,N1)
  170         FORMAT(1X,I19,I7,9I5)
          ELSE
              WRITE (KW,180) (MA(J),J=1,N1)
  180         FORMAT(1X,I19,I7,9I5/(22X,10I5))
          ENDIF
          RETURN
      ENDIF

      IF (JBASE.LE.100000) THEN
          IF (ND.LE.9) THEN
              WRITE (KW,190) (MA(J),J=1,N1)
  190         FORMAT(1X,I19,I8,8I6)
          ELSE
              WRITE (KW,200) (MA(J),J=1,N1)
  200         FORMAT(1X,I19,I8,8I6/(22X,9I6))
          ENDIF
          RETURN
      ENDIF

      IF (ND.LE.4) THEN
          WRITE (KW,210) (MA(J),J=1,N1)
  210     FORMAT(1X,I19,I13,3I11)
      ELSE
          WRITE (KW,220) (MA(J),J=1,N1)
  220     FORMAT(1X,I19,I13,3I11/(22X,4I11))
      ENDIF

END SUBROUTINE

SUBROUTINE FMNTRR(NTR,X,KNAM)

!  Internal routine for trace output of real variables.

!  NTR - 1 for output values
!        2 for input values

!  X   - Double precision value to be printed if NX.EQ.1

!  KNAM - Positive if the routine name is to be printed.

      CHARACTER *4 NAME
      DOUBLE PRECISION X

      IF (NTRACE.EQ.0) RETURN
      IF (NCALL.GT.LVLTRC) RETURN
      IF (NTR.EQ.2 .AND. ABS(NTRACE).EQ.1) RETURN

      IF (NTR.EQ.2 .AND. KNAM.GT.0) THEN
          CALL FMNAME(NAME)
          WRITE (KW,110) NAME
  110     FORMAT(' INPUT TO FM',A4)
      ENDIF
      IF (NTR.EQ.1 .AND. KNAM.GT.0) THEN
          CALL FMNAME(NAME)
          IF (KFLAG.EQ.0) THEN
              WRITE (KW,120) NAME,NCALL,JBASE,NDIG
  120         FORMAT(' FM',A4,15X,'CALL LEVEL =',I2,5X,'JBASE =',  &
                     I10,5X,'NDIG =',I6)
          ELSE
              WRITE (KW,130) NAME,NCALL,JBASE,NDIG,KFLAG
  130         FORMAT(' FM',A4,6X,'CALL LEVEL =',I2,4X,'JBASE =',  &
                     I10,4X,'NDIG =',I6,4X,'KFLAG =',I3)
          ENDIF
      ENDIF

      WRITE (KW,140) X
  140 FORMAT(1X,D30.20)

END SUBROUTINE

SUBROUTINE FMOUT(MA,LINE,LB)

!  Convert a floating multiple precision number to a character array
!  for output.

!  MA   is an FM number to be converted to an A1 character
!       array in base 10 format
!  LINE is the CHARACTER*1 array in which the result is returned.
!  LB   is the length of LINE.

! JFORM1 and JFORM2 (in COMMON) determine the format of LINE.

! JFORM1 = 0  normal setting  ( .314159M+6 )
!        = 1  ES format       ( 3.14159M+5 )
!        = 2  F  format       ( 314159.000 )

! JFORM2 = number of significant digits to display (if JFORM1 = 0, 1)
!        = number of digits after the decimal point (if JFORM1 = 2)

!          If JFORM2.EQ.0 and JFORM1.NE.2 then a default number of
!          digits is chosen.  The default is roughly the full precision
!          of MA.

!          If JFORM2.EQ.0 and JFORM1.EQ.2 then the number is returned in
!          integer format with no decimal point.  Rounding is done as
!          with other settings, so the value displayed is the nearest
!          integer to MA.

!  If JFORM1.EQ.2 and MA is too large or too small to display in the
!  requested format, it is converted using JFORM1=0, JFORM2=0.

!  LINE should be dimensioned at least LOG10(JBASE)*NDIG + 15 on a
!  32-bit machine to allow for up to 10 digit exponents.  Replace
!  15 by 20 if 48-bit integers are used, 25 for 64-bit integers, etc.

      CHARACTER NUMB,NBL,NDPT,M,NPLUS,MINUS,LINE,KSTAR,KCHAR,  &
                NUNKNO(12),NEXPOV(12),NEXPUN(12)

      DIMENSION MA(LUNPCK),LINE(LB),NUMB(10),MD(LUNPCK),MT(LUNPCK),  &
                MS(LUNPCK)

      DATA NUMB /'0','1','2','3','4','5','6','7','8','9'/
      DATA NBL/' '/,NDPT/'.'/,M/'E'/,NPLUS/'+'/,MINUS/'-'/,KSTAR/'*'/
      DATA NUNKNO/' ',' ',' ','U','N','K','N','O','W','N',' ',' '/
      DATA NEXPOV/' ',' ',' ','O','V','E','R','F','L','O','W',' '/
      DATA NEXPUN/' ',' ',' ','U','N','D','E','R','F','L','O','W'/

!             To avoid recursion, FMOUT calls only internal arithmetic
!             routines (FMADD2, FMMPY2, etc.), so no trace printout is
!             done during a call to FMOUT.

      NTRSAV = NTRACE
      NTRACE = 0
      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 32
      DO 110 J = 1, LB
         LINE(J) = NBL
  110 CONTINUE

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO) THEN
          DO 120 J = 1, 12
             LINE(J) = NUNKNO(J)
  120     CONTINUE
          NCALL = NCALL - 1
          NTRACE = NTRSAV
          RETURN
      ENDIF
      IF (MA(1).EQ.KEXPOV) THEN
          DO 130 J = 1, 12
             LINE(J) = NEXPOV(J)
  130     CONTINUE
          LINE(2) = NPLUS
          IF (MA(2).LT.0) LINE(2) = MINUS
          NCALL = NCALL - 1
          NTRACE = NTRSAV
          RETURN
      ENDIF
      IF (MA(1).EQ.KEXPUN) THEN
          DO 140 J = 1, 12
             LINE(J) = NEXPUN(J)
  140     CONTINUE
          LINE(2) = NPLUS
          IF (MA(2).LT.0) LINE(2) = MINUS
          NCALL = NCALL - 1
          NTRACE = NTRSAV
          RETURN
      ENDIF
      KASAVE = KARGSW
      KARGSW = 1
      KRSAVE = KROUND
      KROUND = 1
      JF1SAV = JFORM1
      JF2SAV = JFORM2

!             ND is the number of base 10 digits required.

  150 ND = JFORM2
      IF (JFORM1.EQ.2 .AND. MA(1).GT.0) ND = JFORM2 +  &
          MA(1)*LOG10(REAL(JBASE)) + 1
      IF (ND.LE.1) THEN
          K = NDIG*LOG10(REAL(JBASE))
          ND = MAX(K,JFORM2)
      ENDIF
      IF (JFORM2.LE.0 .AND. JFORM1.LE.1) ND =  &
          1.1 + (NDIG-1)*LOG10(REAL(JBASE))
      IF (ND.LT.2) ND = 2

      IF (LB.LT.ND+6) THEN
          IF (JFORM1.EQ.2) THEN
              JFORM1 = 0
              JFORM2 = 0
              GO TO 150
          ENDIF
          GO TO 370
      ENDIF

!             Convert to the base which is the largest power of 10
!             less than MXBASE and build the output number.

      NPOWER = LOG10(REAL(MXBASE))
      JBSAVE = JBASE
      NDSAVE = NDIG
      MXSAVE = MXEXP
      MXEXP = MXEXP2
      JBASE = 10**NPOWER
      NDIG = ND/NPOWER + 2
      IF (NDIG.LT.2) NDIG = 2
      IF (NDIG.GT.MXNDG2) THEN
          KFLAG = -9
          CALL FMWARN
          GO TO 370
      ENDIF

      IF (MA(2).EQ.0) THEN
          CALL FMIM(0,MS)
          GO TO 210
      ENDIF

!             Check to see if MA is already in a base which is a
!             power of ten.  If so the conversion can be skipped.

      K = NPOWER
      DO 160 J = 1, K
         JBASE = 10**J
         IF (JBASE.EQ.JBSAVE) THEN
             NPOWER = J
             NDIG = ND/NPOWER + 2
             IF (NDIG.LT.2) NDIG = 2
             IF (NDIG.GT.MXNDG2) THEN
                 KFLAG = -9
                 CALL FMWARN
                 GO TO 370
             ENDIF
             CALL FMEQU(MA,MS,NDSAVE,NDIG)
             MS(2) = ABS(MS(2))
             GO TO 210
         ENDIF
  160 CONTINUE

      CALL FMIM(JBSAVE,MD)
      NDS2 = NDSAVE + 1
      CALL FMIM(1,MT)
      KMT = 1

!             Convert the fraction part of MA to the new base.

      KPT = NDS2 + 1
      DO 170 J = 3, NDS2
         KPT = KPT - 1
         IF (MA(KPT).NE.0) GO TO 180
  170 CONTINUE

  180 KEXPSH = KPT - 1
      KDIGIT = ABS(MA(2))
      CALL FMIM(KDIGIT,MS)
      NDIGMS = NDIG

      DO 190 J = 3, KPT
         KDIGIT = MA(J)
         IF (JBSAVE.EQ.2) THEN
             NDIG = MIN(NDIGMS,MAX(2,MS(1)+1))
             CALL FMADD2(MS,MS,MS)
         ELSE
             NDIG = MIN(NDIGMS,MAX(2,MS(1)+MD(1)))
             CALL FMMPY2(MS,MD,MS)
         ENDIF

         IF (KDIGIT.GT.0) THEN
             IF (KMT.NE.KDIGIT) THEN
                 NDIG = MIN(NDIGMS,MAX(2,MD(1)))
                 CALL FMIM(KDIGIT,MT)
                 KMT = KDIGIT
             ENDIF
             NDIG = MIN(NDIGMS,MAX(2,MAX(MS(1),MT(1))+1))
             CALL FMADD2(MS,MT,MS)
         ENDIF
  190 CONTINUE

!             Convert the exponent.

      NDIG = NDIGMS
      CALL FMIM(1,MT)
      K = ABS(MA(1)-KEXPSH)
      IF (MOD(K,2).EQ.1) THEN
          CALL FMEQU(MD,MT,NDIG,NDIG)
      ELSE
          CALL FMIM(1,MT)
      ENDIF

  200 K = K/2
      NDIG = MIN(NDIGMS,MAX(2,MD(1)*2))
      IF (K.GT.0) CALL FMMPY2(MD,MD,MD)
      IF (MOD(K,2).EQ.1) THEN
          NDIG = MIN(NDIGMS,MAX(2,MT(1)+MD(1)))
          CALL FMMPY2(MT,MD,MT)
      ENDIF
      IF (K.GT.1) GO TO 200

      NDIG = NDIGMS
      IF (MA(1)-KEXPSH.LT.0) THEN
          CALL FMDIV2(MS,MT,MS)
      ELSE
          CALL FMMPY2(MS,MT,MS)
      ENDIF

!             Now MS is the value of MA converted to a
!             power of ten base.

!             Next convert it to a character string base 10 for output.

!             KEXP10 is the base 10 exponent.
!             KMS2SD is the number of base 10 significant digits
!                    in MS(2).

  210 MS1 = MS(1)
  220 KEXP10 = NPOWER*MS(1)
      KMS2SD = NPOWER
      K = JBASE
      DO 230 J = 1, NPOWER
         K = K/10
         IF (MS(2).LT.K .AND. MS(2).NE.0) THEN
             KEXP10 = KEXP10 - 1
             KMS2SD = KMS2SD - 1
         ENDIF
  230 CONTINUE

!             For printing using JFORM1 = 1, reduce the exponent to
!             account for the fact that the decimal point and first
!             significant digit will later be swapped.

      IF (JFORM1.EQ.1 .AND. MS(2).NE.0) KEXP10 = KEXP10 - 1

!             Find the position in the unpacked number for rounding.
!             NWORD is the word in which rounding is done, or zero if
!                   no rounding is necessary.
!                   NWORD is set to -1 if JFORM1 is 2 (F format) but no
!                   significant digits would be printed.  This case
!                   defaults to JFORM1 = 0.
!             NVAL gives the position within that word where rounding
!                  occurs.
!             NSD1 is the maximum number of base 10 S.D.'s in NWORD
!                  digits of base 10**NPOWER.
!             NSD2 is the number of base 10 S.D.'s needed to get ND
!                  base 10 digits after the decimal.

      NSD2 = ND
      IF (JFORM1.EQ.2) THEN
          NSD2 = JFORM2 + KEXP10
          IF (NSD2.GT.ND) NSD2 = ND
          NWORD = (NSD2-KMS2SD-1+NPOWER)/NPOWER + 2
          IF (NWORD.LT.2) NWORD = -1
          IF (NWORD.GT.NDIG) NWORD = 0
          IF (NWORD.GE.2 .AND. NSD2.LE.0) NWORD = -1
      ELSE
          NWORD = (ND-KMS2SD-1+NPOWER)/NPOWER + 2
      ENDIF
      NSD1 = KMS2SD + NPOWER*(NWORD-2)
      IF (NWORD.LT.2) THEN
          NVAL = 0
      ELSE
          NVAL = 10**(NSD1-NSD2)
      ENDIF

!             Now do the base 10 rounding.

      IF (NWORD.GE.2) THEN
          X = 0.0
          IF (NVAL.GT.1) X = MOD(MS(NWORD),NVAL)
          IF (NWORD.LT.NDIG+1) X = X + MS(NWORD+1)/REAL(JBASE)
          X = X/NVAL
          IF (X.LT.0.5) GO TO 240
          MS2 = MS(2)
          MS(NWORD) = (MS(NWORD)/NVAL)*NVAL
          MS(NWORD+1) = 0
          MS(NWORD+2) = 0
          MS(NWORD) = MS(NWORD) + NVAL
          IF (MS(NWORD).GE.JBASE) THEN
              NWORD1 = NWORD - 1
              NWORD2 = NWORD - 2
              IF (NWORD.GT.2) THEN
                  CALL FMEQU(MS,MS,NWORD1,NWORD2)
              ELSE
                  MS(1) = MS(1) + 1
                  MS(2) = MS(2)/JBASE
                  MS(3) = 0
              ENDIF
          ENDIF
          IF (MS(1).NE.MS1 .OR. MS(2).NE.MS2) GO TO 220
      ENDIF

!             Build the base 10 character string.

  240 IF (MA(2).LT.0) LINE(1) = MINUS
      LINE(2) = NDPT
      K = 10**KMS2SD
      L = 2
      IF (NWORD.EQ.-1) NSD2 = ND
      DO 250 J = 1, NSD2
         K = K/10
         IF (K.EQ.0) THEN
             K = JBASE/10
             L = L + 1
         ENDIF
         KDIGIT = MS(L)/K
         MS(L) = MOD(MS(L),K)
         LINE(J+2) = NUMB(KDIGIT+1)
  250 CONTINUE

      KA = NSD2 + 3
      KB = ND + 2
      IF (KB.GE.KA) THEN
          DO 260 J = KA, KB
             LINE(J) = NUMB(1)
  260     CONTINUE
      ENDIF

      LINE(ND+3) = M
      LINE(ND+4) = NPLUS
      IF (KEXP10.LT.0) LINE(ND+4) = MINUS
      IF (MA(2).EQ.0) LINE(ND+4) = NBL

!             Build the digits of the base 10 exponent backwards,
!             then reverse them.

      NDE = 1
      KEXP = ABS(KEXP10)
      DO 280 J = 1, LB
         KDIGIT = MOD(KEXP,10)
         LINE(ND+4+J) = NUMB(KDIGIT+1)
         KEXP = KEXP/10
         IF (KEXP.EQ.0) GO TO 290

         IF (ND+5+J.GT.LB) THEN
             DO 270 K = 1, LB
                LINE(K) = KSTAR
  270        CONTINUE
             GO TO 310
         ENDIF

         NDE = NDE + 1
  280 CONTINUE

  290 NDE2 = NDE/2
      IF (NDE2.LT.1) GO TO 310
      K1 = ND + 4
      K2 = ND + 5 + NDE
      DO 300 J = 1, NDE2
         K1 = K1 + 1
         K2 = K2 - 1
         KCHAR = LINE(K1)
         LINE(K1) = LINE(K2)
         LINE(K2) = KCHAR
  300 CONTINUE

!             If JFORM1 is 1 put the first digit left of the decimal.

  310 IF (JFORM1.EQ.1) THEN
          KCHAR = LINE(2)
          LINE(2) = LINE(3)
          LINE(3) = KCHAR
      ENDIF

!             If JFORM1 is 2 put the number into fixed format.

      IF (JFORM1.EQ.2 .AND. JFORM2.GE.0) THEN
          IF (KEXP10.LE.-JFORM2 .OR. KEXP10+2.GT.LB) THEN
              JFORM1 = 0
              JFORM2 = 0
              JBASE = JBSAVE
              NDIG = NDSAVE
              MXEXP = MXSAVE
              DO 320 J = 1, LB
                 LINE(J) = NBL
  320         CONTINUE
              GO TO 150
          ENDIF
          KA = ND + 3
          DO 330 J = KA, LB
             LINE(J) = NUMB(1)
  330     CONTINUE

          IF (KEXP10.GT.0) THEN
              KEXP = KEXP10
              DO 340 J = 1, KEXP
                 LINE(J+1) = LINE(J+2)
  340         CONTINUE
              LINE(KEXP+2) = NDPT
          ENDIF

          IF (KEXP10.LT.0) THEN
              KEXP = -KEXP10
              KA = 3 + KEXP
              KB = LB + 1
              KC = KB - KEXP
              DO 350 J = KA, LB
                 KB = KB - 1
                 KC = KC - 1
                 LINE(KB) = LINE(KC)
                 LINE(KC) = NUMB(1)
  350         CONTINUE
          ENDIF

          JDPT = 0
          DO 360 J = 1, LB
             IF (LINE(J).EQ.NDPT) JDPT = J
             IF (JDPT.GT.0 .AND. J.GT.JDPT+JFORM2) LINE(J) = NBL
  360     CONTINUE
          IF (JFORM2.EQ.0 .AND. JDPT.GT.0) LINE(KEXP+2) = NBL

      ENDIF

!             Restore values and return

      GO TO 390

!             LINE is not big enough to hold the number
!             of digits specified.

  370 KFLAG = -8
      DO 380 J = 1, LB
         LINE(J) = KSTAR
  380 CONTINUE
      CALL FMWARN

  390 JBASE = JBSAVE
      NDIG = NDSAVE
      MXEXP = MXSAVE
      NCALL = NCALL - 1
      KARGSW = KASAVE
      NTRACE = NTRSAV
      KROUND = KRSAVE
      JFORM1 = JF1SAV
      JFORM2 = JF2SAV
END SUBROUTINE

SUBROUTINE FMPI(MA)

!  MA = pi

      DIMENSION MA(LUNPCK)

!             Scratch array usage during FMPI:   M01 - M04

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 33
      IF (ABS(NTRACE).GE.2 .AND. NCALL.LE.LVLTRC) THEN
          WRITE (KW,110)
  110     FORMAT(' INPUT TO FMPI')
      ENDIF
      KASAVE = KARGSW
      KARGSW = 1

!             Increase the working precision.

      NDSAVE = NDIG
      IF (NCALL.EQ.1) THEN
          B = JBASE
          K = 5.0*LOG(10.0)/LOG(B) + 2.0
          NDIG = MAX(NDIG+K,2)
          IF (NDIG.GT.MXNDG2) THEN
              KFLAG = -9
              CALL FMWARN
              MA(1) = KUNKNO
              MA(2) = 1
              DO 120 J = 2, NDSAVE
                 MA(J+1) = 0
  120         CONTINUE
              GO TO 130
          ENDIF
      ENDIF

!             Check to see if pi has previously been computed
!             in base JBASE with sufficient precision.

      IF (NJBPI.EQ.JBASE .AND. NDIGPI.GE.NDIG) THEN
          CALL FMEQU(MPISAV,MA,NDIGPI,NDSAVE)
      ELSE
          NDSV = NDIG
          NDIG = MIN(NDIG+2,MXNDG2)
          CALL FMPI2(MPISAV)
          NJBPI = JBASE
          NDIGPI = NDIG
          CALL FMEQU(MPISAV,MA,NDIG,NDSAVE)
          NDIG = NDSV
      ENDIF

  130 NDIG = NDSAVE
      KARGSW = KASAVE
      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMPI2(MPI)

!  Internal routine to compute pi.
!  The formula used is due to S. Ramanujan:
!                                                (4n)!(1103+26390n)
!  1/pi = (sqrt(8)/9801) * sum(n=0 to infinity) --------------------
!                                               ((n!)**4)(396**(4n))
!  The result is returned in MPI.

      DOUBLE PRECISION X
      DIMENSION MPI(LUNPCK),NSTACK(19)

!             Scratch array usage during FMPI2:   M01 - M04

      NDSAVE = NDIG
      N = -1
      CALL FMI2M(1103,MPI)
      CALL FMI2M(1,M02)
      CALL FMI2M(26390,M03)
      CALL FMI2M(1103,M04)
      MXB2J = MXBASE**2/JBASE

  110 N = N + 1
      LARGE = MXB2J/(4*N + 3)
      J = 4*N + 1
      IF (J.GT.LARGE) THEN
          CALL FMMPYI(M02,J,M02)
          J = J + 1
          CALL FMMPYI(M02,J,M02)
          J = J + 1
          CALL FMMPYI(M02,J,M02)
      ELSE IF (J*(J+1).GT.LARGE) THEN
          K = J*(J+1)
          CALL FMMPYI(M02,K,M02)
          J = J + 2
          CALL FMMPYI(M02,J,M02)
      ELSE
          K = J*(J+1)*(J+2)
          CALL FMMPYI(M02,K,M02)
      ENDIF

      J = N + 1
      LARGE = MXBASE/J
      IF (J.GT.LARGE) THEN
          CALL FMDIVI(M02,J,M02)
          CALL FMDIVI(M02,J,M02)
          CALL FMDIVI(M02,J,M02)
      ELSE IF (J*J.GT.LARGE) THEN
          K = J*J
          CALL FMDIVI(M02,K,M02)
          CALL FMDIVI(M02,J,M02)
      ELSE
          K = J*J*J
          CALL FMDIVI(M02,K,M02)
      ENDIF

!             Break 4/396**4 into 1/(2178*2178*1296).

      J = 2178
      LARGE = MXBASE/J
      IF (J.GT.LARGE) THEN
          CALL FMDIVI(M02,J,M02)
          CALL FMDIVI(M02,J,M02)
          CALL FMDIVI(M02,1296,M02)
      ELSE
          K = J*J
          CALL FMDIVI(M02,K,M02)
          CALL FMDIVI(M02,1296,M02)
      ENDIF

      CALL FMADD(M03,M04,M04)
      CALL FMMPY(M02,M04,M01)

      NDIG = NDSAVE
      CALL FMADD(MPI,M01,MPI)
      NDIG = MAX(2,NDSAVE - (MPI(1) - M01(1)))
      IF (KFLAG.NE.1) GO TO 110
      NDIG = NDSAVE

      CALL FMI2M(8,M02)
      X = 8.0
      X = SQRT(X)
      CALL FMDP2M(X,M04)
      CALL FMDIG(NSTACK,KST)
      DO 120 J = 1, KST
         NDIG = NSTACK(J)
         CALL FMDIV(M02,M04,M01)
         CALL FMADD(M04,M01,M04)
         CALL FMDIVI(M04,2,M04)
  120 CONTINUE
      CALL FMI2M(9801,M03)
      CALL FMMPY(MPI,M04,MPI)
      CALL FMDIV(M03,MPI,MPI)

END SUBROUTINE

SUBROUTINE FMPRNT(MA)

!  Print MA in base 10 format.

!  FMPRNT can be called directly by the user for easy output
!  in M format.  MA is converted using FMOUT and printed.

      DIMENSION MA(LUNPCK)

      KSAVE = KFLAG
      ND = NDIG*LOG10(REAL(JBASE)) + 1
      IF (ND.LT.2) ND = 2
      NEXP = 2.0*LOG10(REAL(MXBASE)) + 6
      LB = MAX(JFORM2+NEXP,ND+NEXP)
      LB = MIN(LB,LMBUFF)
      CALL FMOUT(MA,MBUFF,LB)
      KFLAG = KSAVE
      LAST = LB + 1
      DO 120 J = 1, LB
         IF (MBUFF(LAST-J).NE.' ' .OR. J.EQ.LB) THEN
             L = LAST - J
             WRITE (KW,110) (MBUFF(K),K=1,L)
  110        FORMAT (TR6,73A1)
!  110        FORMAT (6X,73A1)
             RETURN
         ENDIF
  120 CONTINUE
END SUBROUTINE

SUBROUTINE FMPWR(MA,MB,MC)

!  MC = MA ** MB

!  If MB can be expressed exactly as a one word integer, then FMIPWR is
!  used.  This is much faster when MB is small, and using FMIPWR allows
!  MA to be negative.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

!             Scratch array usage during FMPWR:   M01 - M06

!             Convert MB to an integer before changing NDIG.

      KWSV = KWARN
      KWARN = 0
      CALL FMMI(MB,INTMB)
      KWARN = KWSV
      KFL = KFLAG

      CALL FMENTR(34,MA,MB,2,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

!             If the exponent is large, raise the precision.

      IEXTRA = MAX(0,MB(1))
      NDIG = NDIG + IEXTRA
      IF (NDIG.GT.MXNDG2) THEN
          KFLAG = -9
          CALL FMWARN
          MC(1) = KUNKNO
          MC(2) = 1
          DO 110 J = 2, NDSAVE
             MC(J+1) = 0
  110     CONTINUE
          CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
          RETURN
      ENDIF

!             If the exponent is a small integer, call FMIPWR.

      KWSV = KWARN
      KWARN = 0
      CALL FMEQU(MA,M06,NDSAVE,NDIG)
      IF (KFL.EQ.0) THEN
          CALL FMIPWR(M06,INTMB,MC)
      ELSE IF (M06(2).LE.0) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
      ELSE
          CALL FMLN(M06,M06)
          CALL FMEQU(MB,M02,NDSAVE,NDIG)
          CALL FMMPY(M06,M02,M06)
          CALL FMEXP(M06,MC)
      ENDIF
      KWARN = KWSV

!             Round the result and return.

      CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMRDC(MA,MB,JSIN,JCOS,JSWAP)

!  Reduce MA using various trigonometric identities to an equivalent
!  angle MB between 0 and 45 degrees.  The reduction is done in radians
!  if KRAD (in common /FMUSER/) is 1, in degrees if KRAD is 0.
!  JSIN and JCOS are returned +1 or -1 and JSWAP is returned to indicate
!  that the sin and cos functions have been interchanged as follows:

!  JSWAP = 0 means   SIN(MA) = JSIN*SIN(MB)
!                    COS(MA) = JCOS*COS(MB)

!  JSWAP = 1 means   SIN(MA) = JSIN*COS(MB)
!                    COS(MA) = JCOS*SIN(MB)

      DOUBLE PRECISION X
      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMRDC:   M01 - M04

      JSIN = 1
      JCOS = 1
      JSWAP = 0
      NDSAVE = NDIG
      NDIG = NDIG + MAX(0,MA(1))

!             If the argument is too big, return UNKNOWN.

      IF (NDIG.GT.MXNDG2) THEN
          KFLAG = -9
          CALL FMWARN
          MB(1) = KUNKNO
          MB(2) = 1
          DO 110 J = 2, NDSAVE
             MB(J+1) = 0
  110     CONTINUE
          NDIG = NDSAVE
          RETURN
      ENDIF

!             If MA is less than 1/JBASE no reduction is needed.

      IF (MA(1).LT.0) THEN
          NDIG = NDSAVE
          CALL FMEQU(MA,MB,NDIG,NDIG)
          IF (MB(2).LT.0) THEN
              MB(2) = -MB(2)
              JSIN = -1
          ENDIF
          RETURN
      ENDIF

      J = 1
      IF (KRAD.EQ.1) THEN
  120     IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMPI2(MPISAV)
              NJBPI = JBASE
              NDIGPI = NDIG
              NDIG = NDSV
          ENDIF
          CALL FMEQU(MA,M04,NDSAVE,NDIG)
          IF (MA(2).LT.0) JSIN = -1
          M04(2) = ABS(M04(2))
          IF (M04(1).EQ.0) THEN
              CALL FMM2DP(M04,X)
              IF (X.LE.0.75) THEN
                  NDIG = NDSAVE
                  CALL FMEQU(M04,MB,NDIG,NDIG)
                  RETURN
              ENDIF
          ENDIF
          CALL FMMPYI(MPISAV,2,M02)
          IF (FMCOMP(M04,'GE',M02)) THEN
              CALL FMDIV(M04,M02,M01)
              CALL FMINT(M01,M01)
              CALL FMMPY(M01,M02,M01)
              CALL FMSUB(M04,M01,M04)
          ENDIF
          CALL FMEQU(MPISAV,M03,NDIG,NDIG)
          IF (FMCOMP(M04,'GE',M03)) THEN
              JSIN = -JSIN
              CALL FMSUB(M02,M04,M04)
          ENDIF
          CALL FMDIVI(M02,4,M02)
          IF (FMCOMP(M04,'GE',M02)) THEN
              JCOS = -JCOS
              CALL FMSUB(M03,M04,M04)
          ENDIF
          CALL FMDIVI(M03,4,M03)
          IF (FMCOMP(M04,'GE',M03)) THEN
              JSWAP = 1
              CALL FMSUB(M02,M04,M04)
          ENDIF

!             If the reduced argument is close to zero, then
!             cancellation has produced an inaccurate value.
!             Raise NDIG and do the reduction again.

          IF (J.EQ.1 .AND. M04(1).LT.0) THEN
              J = 2
              NDIG = NDIG - M04(1)
              IF (NDIG.GT.MXNDG2) THEN
                  KFLAG = -9
                  CALL FMWARN
                  MB(1) = KUNKNO
                  MB(2) = 1
                  DO 130 J = 2, NDSAVE
                     MB(J+1) = 0
  130             CONTINUE
                  NDIG = NDSAVE
                  RETURN
              ENDIF
              JSIN = 1
              JCOS = 1
              JSWAP = 0
              GO TO 120
          ENDIF

      ELSE

          CALL FMEQU(MA,M04,NDSAVE,NDIG)
          IF (MA(2).LT.0) JSIN = -1
          M04(2) = ABS(M04(2))
          IF (M04(1).EQ.0) THEN
              CALL FMM2DP(M04,X)
              IF (X.LE.44.0) THEN
                  NDIG = NDSAVE
                  CALL FMEQU(M04,MB,NDIG,NDIG)
                  RETURN
              ENDIF
          ENDIF
          CALL FMI2M(360,M02)
          IF (FMCOMP(M04,'GE',M02)) THEN
              CALL FMDIV(M04,M02,M01)
              CALL FMINT(M01,M01)
              CALL FMMPY(M01,M02,M01)
              CALL FMSUB(M04,M01,M04)
          ENDIF
          CALL FMI2M(180,M03)
          IF (FMCOMP(M04,'GE',M03)) THEN
              JSIN = -JSIN
              CALL FMSUB(M02,M04,M04)
          ENDIF
          CALL FMI2M(90,M02)
          IF (FMCOMP(M04,'GE',M02)) THEN
              JCOS = -JCOS
              CALL FMSUB(M03,M04,M04)
          ENDIF
          CALL FMI2M(45,M03)
          IF (FMCOMP(M04,'GE',M03)) THEN
              JSWAP = 1
              CALL FMSUB(M02,M04,M04)
          ENDIF

      ENDIF

!             Round the result and return.

      CALL FMEQU(M04,MB,NDIG,NDSAVE)
      NDIG = NDSAVE
END SUBROUTINE

SUBROUTINE FMRND(ND,NGUARD,KSHIFT)

!  Round MWA to ND digits (base JBASE).

!  MWA is non-negative and has ND+NGUARD+KSHIFT digits.

!  NGUARD is the number of guard digits carried.
!  KSHIFT is 1 if a left shift is pending when MWA(2)=0.

!  Round to position MWA(ND+1+KSHIFT) using the guard digits
!  MWA(ND+2+KSHIFT), ..., MWA(ND+1+NGUARD+KSHIFT).

!  This routine is designed to be called only from within the FM
!  package.  The user should call FMEQU to round numbers.


      IF (KROUND.EQ.0 .AND. NCALL.LE.1) RETURN
      L = ND + 2 + KSHIFT
      IF (2*(MWA(L)+1).LT.JBASE) RETURN
      IF (2*MWA(L).GT.JBASE) THEN
          MWA(L-1) = MWA(L-1) + 1
          MWA(L) = 0
          IF (MWA(L-1).LT.JBASE) RETURN
          GO TO 140
      ENDIF

!             If the first guard digit gives a value close to 1/2 then
!             further guard digits must be examined.

      IF (MOD(JBASE,2).EQ.0) THEN
          IF (2*MWA(L).LT.JBASE) RETURN
          IF (2*MWA(L).EQ.JBASE) THEN
              IF (NGUARD.GE.2) THEN
                  DO 110 J = 2, NGUARD
                     IF (MWA(L+J-1).GT.0) GO TO 130
  110             CONTINUE
              ENDIF

!                       Round to even.

              IF (MOD(MWA(L-1),2).EQ.0) RETURN
          ENDIF
      ELSE
          IF (2*MWA(L)+1.EQ.JBASE) THEN
              IF (NGUARD.GE.2) THEN
                  DO 120 J = 2, NGUARD
                     IF (2*(MWA(L+J-1)+1).LT.JBASE) RETURN
                     IF (2*MWA(L+J-1).GT.JBASE) GO TO 130
  120             CONTINUE
                  RETURN
              ENDIF
          ENDIF
      ENDIF

  130 MWA(L-1) = MWA(L-1) + 1
      MWA(L) = 0

!             Check whether there was a carry in the rounded digit.

  140 KB = L - 1
      IF (KB.GE.3) THEN
          K = KB + 1
          DO 150 J = 3, KB
             K = K - 1
             IF (MWA(K).LT.JBASE) RETURN
             MWA(K-1) = MWA(K-1) + MWA(K)/JBASE
             MWA(K) = MOD(MWA(K),JBASE)
  150     CONTINUE
      ENDIF

!             If there is a carry in the first digit then the exponent
!             must be adjusted and the number shifted right.

      IF (MWA(2).GE.JBASE) THEN
          IF (KB.GE.4) THEN
              K = KB + 1
              DO 160 J = 4, KB
                 K = K - 1
                 MWA(K) = MWA(K-1)
  160         CONTINUE
          ENDIF

          IF (KB.GE.3) MWA(3) = MOD(MWA(2),JBASE)
          MWA(2) = MWA(2)/JBASE
          MWA(1) = MWA(1) + 1
      ENDIF

      IF (KSHIFT.EQ.1 .AND. MWA(2).NE.0) KSHIFT = 0

END SUBROUTINE

SUBROUTINE FMRSLT(MA,MB,MC,KRESLT)

!  Handle results which are special cases, such as overflow,
!  underflow, unknown, etc.

!  MA and MB are the input arguments to an FM subroutine

!  MC is the result which is returned

!  KRESLT is the result code from FMARGS.  Result codes handled here:

!   0 - Perform the normal operation
!   1 - The result is the first input argument
!   2 - The result is the second input argument
!   3 - The result is -OVERFLOW
!   4 - The result is +OVERFLOW
!   5 - The result is -UNDERFLOW
!   6 - The result is +UNDERFLOW
!   7 - The result is -1.0
!   8 - The result is +1.0
!  11 - The result is 0.0
!  12 - The result is 'UNKNOWN'

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      KFSAVE = KFLAG
      IF (KRESLT.EQ.1) THEN
          CALL FMEQU(MA,MC,NDIG,NDIG)
          IF (KSTACK(NCALL).EQ.3 .OR. KSTACK(NCALL).EQ.40) THEN
              KFLAG = 1
          ELSE
              KFLAG = KFSAVE
          ENDIF
          RETURN
      ENDIF

      IF (KRESLT.EQ.2) THEN
          CALL FMEQU(MB,MC,NDIG,NDIG)
          IF (KSTACK(NCALL).EQ.3 .OR. KSTACK(NCALL).EQ.40) THEN
              KFLAG = 1
          ELSE
              KFLAG = KFSAVE
          ENDIF
          RETURN
      ENDIF

      IF (KRESLT.EQ.3 .OR. KRESLT.EQ.4) THEN
          CALL FMIM(0,MC)
          MC(1) = KEXPOV
          MC(2) = 1
          IF (KRESLT.EQ.3) MC(2) = -1
          KFLAG = KFSAVE
          RETURN
      ENDIF

      IF (KRESLT.EQ.5 .OR. KRESLT.EQ.6) THEN
          CALL FMIM(0,MC)
          MC(1) = KEXPUN
          MC(2) = 1
          IF (KRESLT.EQ.5) MC(2) = -1
          KFLAG = KFSAVE
          RETURN
      ENDIF

      IF (KRESLT.EQ.7) THEN
          CALL FMIM(-1,MC)
          KFLAG = KFSAVE
          RETURN
      ENDIF

      IF (KRESLT.EQ.8) THEN
          CALL FMIM(1,MC)
          KFLAG = KFSAVE
          RETURN
      ENDIF

      IF (KRESLT.EQ.11) THEN
          CALL FMIM(0,MC)
          KFLAG = KFSAVE
          RETURN
      ENDIF

      IF (KRESLT.EQ.12 .OR. KRESLT.LT.0 .OR. KRESLT.GT.15) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = KFSAVE
          RETURN
      ENDIF

END SUBROUTINE

SUBROUTINE FMSIGN(MA,MB,MC)

!  MC = SIGN(MA,MB)

!  MC is set to ABS(MA) if MB is positive or zero,
!     or -ABS(MA) if MB is negative.

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 35
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      KWSV = KWARN
      KWARN = 0
      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
      ELSE IF (MB(2).GE.0) THEN
          CALL FMEQU(MA,MC,NDIG,NDIG)
          MC(2) = ABS(MC(2))
      ELSE
          CALL FMEQU(MA,MC,NDIG,NDIG)
          MC(2) = -ABS(MC(2))
      ENDIF

      KWARN = KWSV
      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMSIN(MA,MB)

!  MB = SIN(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMSIN:   M01 - M04

      CALL FMENTR(36,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Reduce the argument, convert to radians if the input is
!             in degrees, and evaluate the function.

      CALL FMRDC(MB,MB,JSIN,JCOS,JSWAP)
      IF (MB(1).EQ.KUNKNO) GO TO 110
      IF (KRAD.EQ.0) THEN
          IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMPI2(MPISAV)
              NJBPI = JBASE
              NDIGPI = NDIG
              NDIG = NDSV
          ENDIF
          CALL FMMPY(MB,MPISAV,MB)
          CALL FMDIVI(MB,180,MB)
      ENDIF
      IF (MB(1).NE.KUNKNO) THEN
          IF (JSWAP.EQ.0) THEN
              CALL FMSIN2(MB,MB)
          ELSE
              CALL FMCOS2(MB,MB)
          ENDIF
      ENDIF

!             Append the sign, round, and return.

      IF (JSIN.EQ.-1 .AND. MB(1).NE.KUNKNO) MB(2) = -MB(2)
  110 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMSIN2(MA,MB)

!  Internal subroutine for MB = SIN(MA) where 0.LE.MA.LE.1.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMSIN2:   M01 - M03

!             LJSUMS = 8*LUNPCK allows for up to eight concurrent sums.
!             Increasing this value will begin to improve the speed of
!             SIN when the base is large and precision exceeds about
!             1,500 decimal digits.

      IF (MA(2).EQ.0) THEN
          CALL FMI2M(0,MA)
          RETURN
      ENDIF
      NDSAVE = NDIG
      KWSV = KWARN
      KWARN = 0

!             Use the direct series
!                  SIN(X) = X - X**3/3! + X**5/5! - ...

!             The argument will be divided by 3**K2 before the series
!             is summed.  The series will be added as J2 concurrent
!             series.  The approximately optimal values of K2 and J2
!             are now computed to try to minimize the time required.
!             N2/2 is the approximate number of terms of the series
!             which will be needed, and L2 guard digits will be carried.

      B = JBASE
      K = 5.0*LOG(10.0)/LOG(B) + 2.0
      T = MAX(NDIG-K,2)
      ALOGB = LOG(B)
      ALOG3 = LOG(3.0)
      ALOGT = LOG(T)
      J2 = 0.081*ALOGB*T**0.3333 + 1.85
      J2 = MAX(1,MIN(J2,LJSUMS/MXNDG2))
      K2 = 0.1*SQRT(T*ALOGB/REAL(J2)) - 0.09*ALOGT + 0.5

      L = -(REAL(MA(1))*ALOGB+LOG(REAL(MA(2))/B +  &
            REAL(MA(3))/(B*B)))/ALOG3 - 0.3
      K2 = K2 - L
      IF (K2.LT.0) THEN
          K2 = 0
          J2 = .43*SQRT(T*ALOGB/(ALOGT+REAL(L)*ALOG3)) + .33
      ENDIF
      IF (J2.LE.1) J2 = 1

      N2 = T*ALOGB/(ALOGT+REAL(L)*ALOG3)
      L2 = LOG(REAL(N2)+3.0**K2)/ALOGB
      NDIG = NDIG + L2
      IF (NDIG.GT.MXNDG2) THEN
          KFLAG = -9
          CALL FMWARN
          MB(1) = KUNKNO
          MB(2) = 1
          DO 110 J = 2, NDSAVE
             MB(J+1) = 0
  110     CONTINUE
          NDIG = NDSAVE
          KWARN = KWSV
          RETURN
      ENDIF
      NDSAV1 = NDIG

!             Divide the argument by 3**K2.

      CALL FMEQU(MA,M02,NDSAVE,NDIG)
      KTHREE = 1
      MAXVAL = MXBASE/3
      IF (K2.GT.0) THEN
          DO 120 J = 1, K2
             KTHREE = 3*KTHREE
             IF (KTHREE.GT.MAXVAL) THEN
                 CALL FMDIVI(M02,KTHREE,M02)
                 KTHREE = 1
             ENDIF
  120     CONTINUE
          IF (KTHREE.GT.1) CALL FMDIVI(M02,KTHREE,M02)
      ENDIF

!             Split into J2 concurrent sums and reduce NDIG while
!             computing each term in the sum as the terms get smaller.

      CALL FMEQU(M02,M03,NDIG,NDIG)
      NTERM = 1
      DO 130 J = 1, J2
         NBOT = NTERM*(NTERM-1)
         IF (NBOT.GT.1) CALL FMDIVI(M03,NBOT,M03)
         NTERM = NTERM + 2
         KPT = (J-1)*(NDIG+1) + 1
         CALL FMEQU(M03,MJSUMS(KPT),NDIG,NDIG)
         M03(2) = -M03(2)
  130 CONTINUE
      CALL FMMPY(M02,M02,M02)
      CALL FMIPWR(M02,J2,MB)

  140 CALL FMMPY(M03,MB,M03)
      DO 150 J = 1, J2
         NBOT = NTERM*(NTERM-1)
         IF (NTERM.GT.MXBASE .OR. NBOT.GT.MXBASE) THEN
             CALL FMDIVI(M03,NTERM,M03)
             NBOT = NTERM - 1
             CALL FMDIVI(M03,NBOT,M03)
         ELSE
             CALL FMDIVI(M03,NBOT,M03)
         ENDIF
         KPT = (J-1)*(NDSAV1+1) + 1
         NDIG = NDSAV1
         CALL FMADD(MJSUMS(KPT),M03,MJSUMS(KPT))
         IF (KFLAG.EQ.1) GO TO 160
         NDIG = NDSAV1 - (MJSUMS(KPT)-M03(1))
         IF (NDIG.LT.2) NDIG = 2
         M03(2) = -M03(2)
         NTERM = NTERM + 2
  150 CONTINUE
      GO TO 140

!             Next put the J2 separate sums back together.

  160 KFLAG = 0
      KPT = (J2-1)*(NDIG+1) + 1
      CALL FMEQU(MJSUMS(KPT),MB,NDIG,NDIG)
      IF (J2.GE.2) THEN
          DO 170 J = 2, J2
             CALL FMMPY(M02,MB,MB)
             KPT = (J2-J)*(NDIG+1) + 1
             CALL FMADD(MB,MJSUMS(KPT),MB)
  170     CONTINUE
      ENDIF

!             Reverse the effect of reducing the argument to
!             compute SIN(MA).

      NDIG = NDSAV1
      IF (K2.GT.0) THEN
          CALL FMI2M(3,M02)
          DO 180 J = 1, K2
             CALL FMMPY(MB,MB,M03)
             CALL FMMPYI(M03,-4,M03)
             CALL FMADD(M02,M03,M03)
             CALL FMMPY(M03,MB,MB)
  180     CONTINUE
      ENDIF

      CALL FMEQU(MB,MB,NDSAV1,NDSAVE)
      NDIG = NDSAVE
      KWARN = KWSV

END SUBROUTINE

SUBROUTINE FMSINH(MA,MB)

!  MB = SINH(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMSINH:   M01 - M03

      CALL FMENTR(37,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      MA1 = MA(1)
      MA2 = MA(2)
      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Work with a positive argument.

      MB(2) = ABS(MB(2))
      KARGSW = 0
      KWRNSV = KWARN
      KWARN = 0

      CALL FMEXP(MB,MB)
      IF (MB(1).EQ.KEXPOV) THEN
          GO TO 110
      ENDIF
      IF (MB(1).EQ.KEXPUN) THEN
          MB(1) = KEXPOV
          MB(2) = -1
          GO TO 110
      ENDIF
      CALL FMI2M(1,M01)

!             If MA is small use the accurate value of EXP(MA) - 1
!             which is computed by FMEXP and is stored in M03 to
!             avoid cancellation error.

      IF (MA1.GT.0) THEN
          CALL FMDIV(M01,MB,M01)
          CALL FMSUB(MB,M01,MB)
      ELSE
          CALL FMADD(MB,M01,M01)
          CALL FMMPY(M03,M01,M01)
          CALL FMDIV(M01,MB,MB)
      ENDIF
      CALL FMDIVI(MB,2,MB)

!             Round and return.

  110 IF (MA2.LT.0) MB(2) = -MB(2)
      KWARN = KWRNSV
      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMSP2M(X,MA)

!  MA = X

!  Convert a single precision number to FM format.

!  In general the relative accuracy of the number returned is only
!  the relative accuracy of a machine precision number.  This may be
!  true even if X can be represented exactly in the machine floating
!  point number system.

      DOUBLE PRECISION Y
      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 38
      Y = X
      IF (NTRACE.NE.0) CALL FMNTRR(2,Y,1)

      CALL FMDM(Y,MA)

      IF (NTRACE.NE.0) CALL FMNTR(1,MA,MA,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMSQRT(MA,MB)

!  MB = SQRT(MA)

      DOUBLE PRECISION X,XB
      DIMENSION MA(LUNPCK),MB(LUNPCK),NSTACK(19)

!             Scratch array usage during FMSQRT:   M01 - M02

      CALL FMENTR(39,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      MA1 = MA(1)
      CALL FMEQU(MA,M02,NDSAVE,NDIG)

!             Generate the first approximation.

      M02(1) = 0
      CALL FMM2DP(M02,X)
      X = SQRT(X)
      KE = MA1/2
      IF (MOD(ABS(MA1),2).EQ.1) THEN
          XB = JBASE
          X = X*SQRT(XB)
          KE = (MA1-1)/2
      ENDIF
      CALL FMDP2M(X,MB)
      MB(1) = MB(1) + KE

!             Initialize.

      M02(1) = MA1
      CALL FMDIG(NSTACK,KST)

!             Newton iteration.

      DO 110 J = 1, KST
         NDIG = NSTACK(J)
         CALL FMDIV(M02,MB,M01)
         CALL FMADD(MB,M01,MB)
         CALL FMDIVI(MB,2,MB)
  110 CONTINUE

!             Round the result and return.

      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMSUB(MA,MB,MC)

!  MC = MA - MB

      DIMENSION MA(LUNPCK),MB(LUNPCK),MC(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 40
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MB,2)

      KFLG1 = 0
      IF (MB(1).GT.MA(1) .OR. MA(2).EQ.0) KFLG1 = 1
      IF (MB(2).EQ.0) KFLG1 = 0
      K1 = MB(2)
      MB(2) = -MB(2)

      CALL FMADD2(MA,MB,MC)

!             If MB and MC are distinct arrays then the sign of MB must
!             be restored.  If they refer to identical memory locations
!             then MB now contains the result of the subtraction, so it
!             should not be restored.  Detecting which of these two
!             cases is true is tricky because some optimizing compilers
!             remove lines 1 and 3 of the following:
!                   K2 = MC(2)
!                   MB(2) = K1
!                   MC(2) = K2
!             The code used below tries to do the same thing in an
!             obscure way which defeats the optimizer.
!             If this still fails for some compiler, copy MB to a local
!             array and negate that as the argument to FMADD2.  This
!             adds one unpacked array to the program size and makes
!             subtraction run noticeably slower than addition.
!             See the Portability Notes before FMSET for a test to
!             detect this failure.

      IF (MB(2).NE.MC(2)) THEN
          MB(2) = K1
      ELSE
          K2 = (MB(2)+MC(2))/2
          MB(2) = K1
          MC(2) = K2
      ENDIF

!             If MA was smaller than MB then KFLAG = 1 returned from
!             FMADD means the result from FMSUB is the opposite of the
!             input argument of larger magnitude, so reset KFLAG.

      IF (KFLAG.EQ.1 .AND. KFLG1.EQ.1) KFLAG = 0

      IF (NTRACE.NE.0) CALL FMNTR(1,MC,MC,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMTAN(MA,MB)

!  MB = TAN(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMTAN:   M01 - M04

      CALL FMENTR(41,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Reduce the argument, convert to radians if the input is
!             in degrees, and evaluate the function.

      CALL FMRDC(MB,MB,JSIN,JCOS,JSWAP)
      IF (MB(1).EQ.KUNKNO) GO TO 110
      IF (MB(2).EQ.0) THEN
          IF (JSWAP.EQ.1) THEN
              CALL FMIM(0,MB)
              MB(1) = KUNKNO
              MB(2) = 1
              KFLAG = -4
              CALL FMWARN
          ENDIF
          GO TO 110
      ENDIF
      IF (KRAD.EQ.0) THEN
          IF (NJBPI.NE.JBASE .OR. NDIGPI.LT.NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,MXNDG2)
              CALL FMPI2(MPISAV)
              NJBPI = JBASE
              NDIGPI = NDIG
              NDIG = NDSV
          ENDIF
          CALL FMMPY(MB,MPISAV,MB)
          CALL FMDIVI(MB,180,MB)
      ENDIF
      IF (MB(1).NE.KUNKNO) THEN
          CALL FMSIN2(MB,MB)
          CALL FMMPY(MB,MB,M03)
          CALL FMI2M(1,M02)
          CALL FMSUB(M02,M03,M03)
          CALL FMSQRT(M03,M03)
          IF (JSWAP.EQ.0) THEN
              CALL FMDIV(MB,M03,MB)
          ELSE
              CALL FMDIV(M03,MB,MB)
          ENDIF
      ENDIF

!             Append the sign, round, and return.

      IF (JSIN*JCOS.EQ.-1) MB(2) = -MB(2)
  110 CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMTANH(MA,MB)

!  MB = TANH(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

!             Scratch array usage during FMTANH:   M01 - M03

      CALL FMENTR(42,MA,MA,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KRESLT.NE.0) RETURN

      KARGSW = 0
      KWRNSV = KWARN
      KWARN = 0
      MA1 = MA(1)
      MA2 = MA(2)
      CALL FMEQU(MA,MB,NDSAVE,NDIG)

!             Work with a positive argument.

      MB(2) = ABS(MB(2))

      CALL FMEXP(MB,MB)
      IF (MB(1).EQ.KEXPOV) THEN
          CALL FMI2M(1,MB)
          GO TO 110
      ENDIF
      IF (MB(1).EQ.KEXPUN) THEN
          CALL FMI2M(-1,MB)
          GO TO 110
      ENDIF
      CALL FMI2M(1,M01)

!             If MA is small use the accurate value of EXP(MA) - 1
!             which is computed by FMEXP and is stored in M03 to
!             avoid cancellation error.

      IF (MA1.GT.0) THEN
          CALL FMDIV(M01,MB,M01)
          CALL FMSUB(MB,M01,M02)
          CALL FMADD(MB,M01,M01)
          CALL FMDIV(M02,M01,MB)
      ELSE
          CALL FMADD(MB,M01,M01)
          CALL FMMPY(M03,M01,M01)
          CALL FMI2M(2,M02)
          CALL FMADD(M01,M02,MB)
          CALL FMDIV(M01,MB,MB)
      ENDIF

!             Round and return.

  110 IF (MA2.LT.0) MB(2) = -MB(2)
      KWARN = KWRNSV
      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

SUBROUTINE FMTRAP(MA)

!  If MA has overflowed or underflowed, replace it by the appropriate
!  symbol.

      DIMENSION MA(LUNPCK)

      IF (NCALL.LE.0) RETURN
      IF (MA(1).GT.MXEXP+1) THEN
          MA(1) = KEXPOV
          IF (MA(2).GT.0) THEN
              MA(2) = 1
          ELSE
              MA(2) = -1
          ENDIF
          KFLAG = -5
          CALL FMWARN
      ENDIF
      IF (MA(1).LT.-MXEXP) THEN
          MA(1) = KEXPUN
          IF (MA(2).GT.0) THEN
              MA(2) = 1
          ELSE
              MA(2) = -1
          ENDIF
          KFLAG = -6
          CALL FMWARN
      ENDIF

END SUBROUTINE

SUBROUTINE FMULP(MA,MB)

!  MB = The value of one Unit in the Last Place of MA at the current
!       base and precision.

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      KFLAG = 0
      NCALL = NCALL + 1
      KSTACK(NCALL) = 43
      IF (NTRACE.NE.0) CALL FMNTR(2,MA,MA,1)

      MA1 = MA(1)
      N1 = NDIG + 1
      DO 110 J = 3, N1
         MWA(J) = 0
  110 CONTINUE
      MWA(2) = 1
      IF (MA(2).LT.0) MWA(2) = -1
      MWA(1) = MA(1) - NDIG + 1
      IF (MA(2).EQ.0 .OR. MA(1).GE.KEXPOV) THEN
          CALL FMIM(0,MB)
          MB(1) = KUNKNO
          MB(2) = 1
          KFLAG = -4
          IF (MA1.NE.KUNKNO) CALL FMWARN
      ELSE
          KWRNSV = KWARN
          IF (MA1.EQ.KEXPUN) KWARN = 0
          CALL FMMOVE(MB)
          KWARN = KWRNSV
      ENDIF

      IF (NTRACE.NE.0) CALL FMNTR(1,MB,MB,1)
      NCALL = NCALL - 1
END SUBROUTINE

SUBROUTINE FMWARN

!  Called by one of the FM routines to print a warning message
!  if any error condition arises in that routine.

      CHARACTER *4 NAME

      IF (KFLAG.GE.0 .OR. NCALL.NE.1 .OR. KWARN.LE.0) RETURN
      NCS = NCALL
      CALL FMNAME(NAME)
      WRITE(KW,110) KFLAG,NAME
  110 FORMAT (/' ERROR OF TYPE KFLAG =',I3,  &
              ' IN FM PACKAGE IN ROUTINE FM',A4/)

  120 NCALL = NCALL - 1
      IF (NCALL.GT.0) THEN
          CALL FMNAME(NAME)
          WRITE (KW,130) NAME
  130     FORMAT( ' CALLED FROM FM',A4)
          GO TO 120
      ENDIF

      IF (KFLAG.EQ.-1) THEN
          WRITE (KW,140) MXNDIG
  140     FORMAT(' NDIG MUST BE BETWEEN 2 AND',I10/)
      ELSE IF (KFLAG.EQ.-2) THEN
          WRITE (KW,150) MXBASE
  150     FORMAT(' JBASE MUST BE BETWEEN 2 AND',I10/)
      ELSE IF (KFLAG.EQ.-3) THEN
          WRITE (KW,160)
  160     FORMAT(' AN INPUT ARGUMENT IS NOT A VALID FM NUMBER.',  &
                 '  ITS EXPONENT IS OUT OF RANGE.'/)
          WRITE (KW,170)
  170     FORMAT(' UNKNOWN HAS BEEN RETURNED.'/)
      ELSE IF (KFLAG.EQ.-4 .OR. KFLAG.EQ.-7) THEN
          WRITE (KW,180)
  180     FORMAT(' INVALID INPUT ARGUMENT FOR THIS ROUTINE.'/)
          WRITE (KW,170)
      ELSE IF (KFLAG.EQ.-5) THEN
          WRITE (KW,190)
  190     FORMAT(' THE RESULT HAS OVERFLOWED.'/)
      ELSE IF (KFLAG.EQ.-6) THEN
          WRITE (KW,200)
  200     FORMAT(' THE RESULT HAS UNDERFLOWED.'/)
      ELSE IF (KFLAG.EQ.-8) THEN
          WRITE (KW,210)
  210     FORMAT(' THE RESULT ARRAY IS NOT BIG ENOUGH TO HOLD THE',  &
                 ' OUTPUT CHARACTER STRING'/' IN THE CURRENT FORMAT.'/  &
                 ' THE RESULT ''***...***'' HAS BEEN RETURNED.'/)
      ELSE IF (KFLAG.EQ.-9) THEN
          WRITE (KW,220)
  220     FORMAT(' PRECISION COULD NOT BE RAISED ENOUGH TO PROVIDE ALL'  &
                 ,' REQUESTED GUARD DIGITS.'/)
          WRITE (KW,230) NDIG,MXNDG2
  230     FORMAT(' REQUESTED NDIG=',I7,'.  MAXIMUM AVAILABLE NDIG=',I7/)
          WRITE (KW,170)
      ELSE IF (KFLAG.EQ.-10) THEN
          IF (KSTACK(NCS).EQ.25) THEN
              WRITE (KW,240)
  240         FORMAT(' AN FM NUMBER WAS TOO SMALL IN MAGNITUDE TO ',  &
                     'CONVERT TO SINGLE PRECISION.'/)
          ELSE
              WRITE (KW,250)
  250         FORMAT(' AN FM NUMBER WAS TOO SMALL IN MAGNITUDE TO ',  &
                     'CONVERT TO DOUBLE PRECISION.'/)
          ENDIF
          WRITE (KW,260)
  260     FORMAT(' ZERO HAS BEEN RETURNED.'/)
      ENDIF

      NCALL = NCS
      IF (KWARN.GE.2) STOP
END SUBROUTINE

END MODULE XP_REAL_MODULE



!program test_xp
!use xp_real_module
!type(xp_real) :: x, y, z
!call xp_set(100)
!x = 1.0
!y = 4.0
!z = y * atan(x)
!jform2 = 100
!call xp_print(z)
!end program test_xp
