!     FM 1.2                    David M. Smith                  8-17-01


!  The routines in this package perform multiple precision arithmetic and
!  functions on three kinds of numbers.
!  FM routines handle floating-point real multiple precision numbers,
!  IM routines handle integer multiple precision numbers, and
!  ZM routines handle floating-point complex multiple precision numbers.


!  1. INITIALIZING THE PACKAGE

!  The variables that contain values to be shared by the different routines are
!  located in module FMVALS in file FMSAVE.f90.  Variables that are described
!  below for controlling various features of the FM package are found in this
!  module.  They are initialized to default values assuming 32-bit integers and
!  64-bit double precision representation of the arrays holding multiple
!  precision numbers.  The base and number of digits to be used are initialized
!  to give slightly more than 50 decimal digits.  Subroutine FMVARS can be used
!  to get a list of these variables and their values.

!  The intent of module FMVALS is to hide the FM internal variables from the
!  user's program, so that no name conflicts can occur.  Subroutine FMSETVAR can
!  be used to change the variables listed below to new values.  It is not always
!  safe to try to change these variables directly by putting USE FMVALS into the
!  calling program and then changing them by hand.  Some of the saved constants
!  depend upon others, so that changing one variable may cause errors if others
!  depending on that one are not also changed.  FMSETVAR automatically updates
!  any others that depend upon the one being changed.

!  Subroutine FMSET also initializes these variables.  It tries to compute the
!  best value for each, and it checks several of the values set in FMVALS to see
!  that they are reasonable for a given machine.  FMSET can also be called to
!  set or change the current precision level for the multiple precision numbers.

!  Calling FMSET is optional in version 1.2 of the FM package.  In previous
!  versions one call was required before any other routine in the package could
!  be used.

!  The routine ZMSET from version 1.1 is no longer needed, and the complex
!  operations are automatically initialized in FMVALS.  It has been left in the
!  package for compatibility with version 1.1.


!  2.  REPRESENTATION OF FM NUMBERS

!  MBASE is the base in which the arithmetic is done.  MBASE must be
!        bigger than one, and less than or equal to the square root of
!        the largest representable integer.  For best efficiency MBASE
!        should be large, but no more than about 1/4 of the square root
!        of the largest representable integer.  Input and output
!        conversions are much faster when MBASE is a power of ten.

!  NDIG  is the number of base MBASE digits that are carried in the
!        multiple precision numbers.  NDIG must be at least two.  The
!        upper limit for NDIG is defined in FMVALS and is restricted
!        only by the amount of memory available.

!  Sometimes it is useful to dynamically vary NDIG during the program.  Routine
!  FMEQU should be used to round numbers to lower precision or zero-pad them to
!  higher precision when changing NDIG.

!  The default value of MBASE is a large power of ten.  FMSET also sets MBASE to
!  a large power of ten.  For an application where another base is used, such as
!  simulating a given machine's base two arithmetic, use subroutine FMSETVAR to
!  change MBASE, so that the other internal values depending on MBASE will be
!  changed accordingly.

!  There are two representations for a floating point multiple precision number.
!  The unpacked representation used by the routines while doing the computations
!  is base MBASE and is stored in NDIG+3 words.  A packed representation is
!  available to store the numbers in the user's program in compressed form.  In
!  this format, the NDIG (base MBASE) digits of the mantissa are packed two per
!  word to conserve storage.  Thus the external, packed form of a number
!  requires (NDIG+1)/2+3 words.

!  This version uses double precision arrays to hold the numbers.  Version 1.0
!  of FM used integer arrays, which are faster on some machines.  The package
!  can be changed to use integer arrays --- see section 11 on EFFICIENCY below.

!  The unpacked format of a floating multiple precision number is as follows.
!  A number MA is kept in an array with MA(1) containing the exponent, and
!  MA(2) through MA(NDIG+1) each containing one digit of the mantissa, expressed
!  in base MBASE.  The array is dimensioned to start at MA(-1), with the
!  sign of the number (+1 or -1) held in MA(-1), and the  approximate number
!  of bits of precision stored in MA(0).  This precision value is intended to
!  be used by FM functions that need to monitor cancellation error in addition
!  and subtraction.  The cancellation monitor code is usually disabled for user
!  calls, and FM functions only check for cancellation when they must.  Tracking
!  cancellation causes most routines to run slower, with addition and
!  subtraction being affected the most.  The exponent is a power of MBASE 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.

!  In both representations 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.

!  For MBASE = 10,000 and NDIG = 4, if array MA holds the number -pi, it would
!  have these representations:
!                   Word 1         2         3         4         5

!         Unpacked:      1         3      1415      9265      3590
!         Packed:        1     31415  92653590

!  In both formats MA(0) would be 42, indicating that the mantissa has about 42
!  bits of precision, and MA(-1) = -1 since the number is negative.

!  Because of normalization in a large base, the equivalent number of base 10
!  significant digits for an FM number may be as small as
!  LOG10(MBASE)*(NDIG-1) + 1.

!  The integer routines use the FM format to represent numbers, without the
!  number of digits (NDIG) being fixed.  Integers in IM format are essentially
!  variable precision, using the minimum number of words to represent each
!  value.

!  For programs using both FM and IM numbers, FM routines should not be called
!  with IM numbers, and IM routines should not be called with FM numbers, since
!  the implied value of NDIG used for an IM number may not match the explicit
!  NDIG expected by an FM routine.  Use the conversion routines IMFM2I and
!  IMI2FM to change between the FM and IM formats.

!  The format for complex FM numbers (called ZM numbers below) is very similar
!  to that for real FM numbers.  Each ZM array holds two FM numbers to represent
!  the real and imaginary parts of a complex number.  Each ZM array is twice as
!  long as a corresponding FM array, with the imaginary part starting at the
!  midpoint of the array.  As with FM, there are packed and unpacked formats for
!  the numbers.


!  3. INPUT/OUTPUT ROUTINES

!  All versions of the input routines perform free-format conversion from
!  characters to FM numbers.

!  a. Conversion to or from a character array

!     FMINP converts from a character*1 array to an FM number.

!     FMOUT converts an FM number to base 10 and formats it for output as an
!           array of type character*1.  The output is left justified in the
!           array, and the format is defined by two variables in module FMVALS,
!           so that a separate format definition does not have to be provided
!           for each output call.

!     JFORM1 and JFORM2 define a default output format.

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

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

!  b. Conversion to or from a character string

!     FMST2M converts from a character string to an FM number.

!     FMFORM converts an FM number to a character string according to a format
!            provided in each call.  The format description is more like that of
!            a Fortran FORMAT statement, and integer or fixed-point output is
!            right justified.

!  c. Direct read or write

!     FMPRNT uses FMOUT to print one FM number.

!     FMFPRT uses FMFORM to print one FM number.

!     FMWRIT writes FM numbers for later input using FMREAD.

!     FMREAD reads FM numbers written by FMWRIT.

!  The values given to JFORM1 and JFORM2 can be used to define a default output
!  format when FMOUT or FMPRNT are called.  The explicit format used in a call
!  to FMFORM or FMFPRT overrides the settings of JFORM1 and JFORM2.

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

!  For multiple precision integers, the corresponding routines IMINP, IMOUT,
!  IMST2M, IMFORM, IMPRNT, IMFPRT, IMWRIT, and IMREAD provide similar input and
!  output conversions.  For output of IM numbers, JFORM1 and JFORM2 are ignored
!  and integer format (JFORM1=2, JFORM2=0) is used.

!  For ZM numbers, the corresponding routines ZMINP, ZMOUT, ZMST2M, ZMFORM,
!  ZMPRNT, ZMFPRT, ZMWRIT, and ZMREAD provide similar input and output
!  conversions.

!  For the output format of ZM numbers, JFORM1 and JFORM2 determine the default
!  format for the individual parts of a complex number as with FM numbers.

!     JFORMZ determines the combined output format of the real and
!            imaginary parts.

!     JFORMZ = 1  normal setting    :    1.23 - 4.56 i
!            = 2  use capital I     :    1.23 - 4.56 I
!            = 3  parenthesis format:  ( 1.23 , -4.56 )

!     JPRNTZ controls whether to print real and imaginary parts
!            on one line whenever possible.

!     JPRNTZ = 1  print both parts as a single string :
!                     1.23456789M+321 - 9.87654321M-123 i
!            = 2  print on separate lines without the 'i' :
!                     1.23456789M+321
!                    -9.87654321M-123

!  For further description of these routines, see sections 9 and 10 below.


!  4. ARITHMETIC TRACING

!  NTRACE and LVLTRC control trace printout from the package.

!  NTRACE =  0   No output except warnings and errors.  (Default)
!         =  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 MBASE 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 MBASE 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.  Default is 1.

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


!  5. ERROR CONDITIONS

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

!  Subroutine FMFLAG is provided to give the user access to the current
!  condition code.  For example, to set the user's local variable LFLAG
!  to FM's internal KFLAG value:      CALL FMFLAG(LFLAG)

!    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 NDIGMX.
!          = -2     MBASE 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 (array) to FMINP was not legal.
!          = -8     The character array was not large enough in an
!                   input or output routine.
!          = -9     Precision could not be raised enough to provide all
!                        requested guard digits.  Increasing the value
!                        of NDIGMX in file FMSAVE.f90 may fix this.
!                        UNKNOWN was returned.
!          = -10    An FM input argument was too small in magnitude to
!                        convert to the machine's single or double
!                        precision in FMM2SP or FMM2DP.  Check that the
!                        definitions of SPMAX and DPMAX in file FMSAVE.f90
!                        are correct for the current machine.
!                        Zero was returned.
!          = -11    Array MBERN is not dimensioned large enough for the
!                        requested number of Bernoulli numbers.
!          = -12    Array MJSUMS is not dimensioned large enough for
!                        the number of coefficients needed in the
!                        reflection formula in FMPGAM.

!  When a negative KFLAG condition is encountered, the value of KWARN
!  determines the action to be taken.

!  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.

!  The default setting is KWARN = 1.

!  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.

!  IM routines handle exceptions like OVERFLOW or UNKNOWN in the same way as FM
!  routines.  When using IMMPY, the product of two large positive integers will
!  return +OVERFLOW.  The routines IMMPYM and IMPMOD can be used to obtain a
!  modular result without overflow.  The largest representable IM integer is
!  MBASE**NDIGMX - 1.  For example, if MBASE is 10**7 and NDIGMX is set to 256,
!  integers less than 10**1792 can be used.


!  6. OTHER OPTIONS

!  KRAD = 0     All angles in the trigonometric functions and inverse functions
!               are measured in degrees.
!       = 1     All angles are measured in radians.  (Default)

!  KROUND = -1  All results are rounded toward minus infinity.
!         =  0  All results are rounded toward zero (chopped).
!         =  1  All results are rounded to the nearest FM number, or to the
!               value with an even last digit if the result is halfway
!               between two FM numbers.  (Default)
!         =  2  All results are rounded toward plus infinity.

!         In all cases, while a function is being computed all intermediate
!         results are rounded to nearest, with only the final result being
!         rounded according to KROUND.

!  KRPERF = 0   A smaller number of guard digits used, to give nearly perfect
!               rounding.  This number is chosen so that the last intermediate
!               result should have error less than 0.001 unit in the last place
!               of the final rounded result.  (Default)
!         = 1   Causes more guard digits to be used, to get perfect rounding in
!               the mode set by KROUND.  This slows execution speed.

!         If a small base is used for the arithmetic, like MBASE = 2, 10, or 16,
!         FM assumes that the arithmetic hardware for some machine is being
!         simulated, so perfect rounding is done without regard for the value
!         of KRPERF.
!         If KROUND = 1, then KRPERF = 1 means returned results are no more than
!         0.500 units in the last place from the exact mathematical result,
!         versus 0.501 for KRPERF = 0.
!         If KROUND is not 1, then KRPERF = 1 means returned results are no more
!         than 1.000 units in the last place from the exact mathematical result,
!         versus 1.001 for KRPERF = 0.

!  KSWIDE defines the maximum screen width to be used for all unit KW output.
!         Default is 80.

!  KESWCH controls the action taken in FMINP and other input routines for
!         strings like 'E7' that have no digits before the exponent field.
!         This is sometimes a convenient abbreviation when doing interactive
!         keyboard input.
!         KESWCH = 1 causes 'E7' to translate like '1.0E+7'.  (Default)
!         KESWCH = 0 causes 'E7' to translate like '0.0E+7' and give 0.

!  CMCHAR defines the exponent letter to be used for FM variable output.
!         Default is 'M', as in 1.2345M+678.
!         Change it to 'E' for output to be read by a non-FM program.

!  KDEBUG = 0   No error checking is done to see if input arguments are valid
!               and parameters like NDIG and MBASE are correct upon entry to
!               each routine.  (Default)
!         = 1   Some error checking is done.  (Slower speed)

!  See module FMVALS in file FMSAVE.f90 for additional description of these and
!  other variables defining various FM conditions.


!  7. ARRAY DIMENSIONS

!  The dimensions of the arrays in the FM package are defined using parameters
!  NDIGMX and NBITS.
!  NDIGMX is the maximum value the user may set for NDIG.
!  NBITS  is the number of bits used to represent integers for a given machine.
!         See the EFFICIENCY discussion below.

!  The standard version of FM sets NDIGMX = 55, so on a 32-bit machine using
!  MBASE = 10**7 the maximum precision is about 7*54+1 = 379 significant
!  digits.  Previous versions of FM set NDIGMX = 256.  Two reasons for making
!  this change are:
!  (a) Almost all applications using FM use only 30 to 50 significant digits
!      for checking double or quadruple precision results, and the larger
!      arrays are wasted space.
!  (b) Most FM applications use the derived type interface so that the number
!      of changes to existing code is minimized.  Many compilers implement the
!      FM interface by doing copy in / copy out argument passing of the derived
!      types.  Copying the entire large array when only a small part of it is
!      being used causes the derived type arithmetic to be slow compared to
!      making direct calls to the subroutines.  Setting NDIGMX to be only
!      slightly higher than a program actually uses minimizes any performance
!      penalty for the derived type arithmetic.

!  To change dimensions so that 10,000 significant digit calculation can be
!  done, NDIGMX needs to be at least  10**4/7 + 5 = 1434.  This allows for a
!  few user guard digits to be defined when the precision is changed using
!  CALL FMSET(10000).  Changing 'NDIGMX = 55' to 'NDIGMX = 1434' in FMSAVE.f90
!  will define all the new array sizes.

!  If NDIG much greater than 256 is to be used and elementary functions will
!  be needed, they will be faster if array MJSUMS is larger.  The parameter
!  defining the size of MJSUMS is set in the standard version by
!      LJSUMS = 8*(LUNPCK+3)
!  The 8 means that up to eight concurrent sums can be used by the elementary
!  functions.  The approximate number needed for best speed is given by
!      0.051*Log(MBASE)*NDIG**0.333 + 1.85
!  For example, with MBASE=10**7 and NDIG=1434 this gives 11.  Changing
!  'LJSUMS = 8*(LUNPCK+3)' to 'LJSUMS = 11*(LUNPCK+3)' in FMSAVE.f90 will give
!  slightly better speed.

!  FM numbers in packed format have dimension -1:LPACK, and those in unpacked
!  format have dimension -1:LUNPCK.

!  The parameters LPACKZ and LUNPKZ define the size of the packed and unpacked
!  ZM arrays.  The real part starts at the beginning of the array, and the
!  imaginary part starts at word KPTIMP for packed format or at word KPTIMU for
!  unpacked format.


!  8. PORTABILITY

!  In FMSET several variables are set to machine-dependent values, and many of
!  the variables initialized in module FMVALS in file FMSAVE.f90 are checked to
!  see that they have reasonable values.  FMSET will print warning messages on
!  unit KW for any of the FMVALS variables that seem to be poorly initialized.

!  If an FM run fails, call FMVARS to get a list of all the FMVALS variables
!  printed on unit KW.  Setting KDEBUG = 1 at the start may also identify some
!  errors.

!  Some compilers object to a function like FMCOMP with side effects such as
!  changing KFLAG or other module variables.  Blocks of code in FMCOMP and
!  IMCOMP that modify these variables are identified so they may be removed or
!  commented out to produce a function without side effects.  This disables
!  trace printing in FMCOMP and IMCOMP, and error codes are not returned in
!  KFLAG.  See FMCOMP and IMCOMP for further details.

!  In FMBER2 and FMPGAM several constants are used that require the machine's
!  integer word size to be at least 32 bits.


!  9.  LIST OF ROUTINES

!  These are the FM routines that are designed to be called by the user.
!  All are subroutines except logical function FMCOMP.
!  MA, MB, MC refer to FM format numbers.

!  In Fortran-90 and later versions of the Fortran standard, it is potentially
!  unsafe to use the same array more than once in the calling sequence.  The
!  operation MA = MA + MB should not be written as
!        CALL FMADD(MA,MB,MA)
!  since the compiler is allowed to pass the three arguments with a copy in /
!  copy out mechanism.  This means the third argument, containing the result,
!  might not be copied out last, and then a later copy out of the original
!  input MA could destroy the computed result.

!  One solution is to use a third array and then put the result back in MA:
!        CALL FMADD(MA,MB,MC)
!        CALL FMEQ(MC,MA)

!  When the first call is doing one of the "fast" operations like addition,
!  the extra call to move the result back to MA can cause a noticeable loss in
!  efficiency.  To avoid this, separate routines are provided for the basic
!  arithmetic operations when the result is to be returned in the same array
!  as one of the inputs.

!  A routine name with a suffix of  "_R1" returns the result in the first
!  input array, and a suffix of "_R2" returns the result in the second input
!  array.  The example above would then be:
!        CALL FMADD_R1(MA,MB)

!  These routines each have one less argument than the original version, since
!  the output is re-directed to one of the inputs.  The result array should
!  not be the same as any input array when the original version of the routine
!  is used.

!  The routines that can be used this way are listed below.  For others, like
!        CALL FMEXP(MA,MA)
!  the relative cost of doing an extra copy is small.  This one should become
!        CALL FMEXP(MA,MB)
!        CALL FMEQ(MB,MA)

!  If the derived-type interface is used, as in
!        TYPE (FM) A,B
!        ...
!        A = A + B
!  there is no problem putting the result back into A, since the interface
!  routine creates a temporary scratch array for the result of A + B, allowing
!  copy in / copy out to work.

!  For each of these routines there is also a version available for which the
!  argument list is the same but all FM numbers are in packed format.  The
!  routines using packed numbers have the same names except 'FM' is replaced by
!  'FP' at the start of each name.


!  FMABS(MA,MB)         MB = ABS(MA)

!  FMACOS(MA,MB)        MB = ACOS(MA)

!  FMADD(MA,MB,MC)      MC = MA + MB

!  FMADD_R1(MA,MB)      MA = MA + MB

!  FMADD_R2(MA,MB)      MB = MA + MB

!  FMADDI(MA,IVAL)      MA = MA + IVAL   Increment an FM number by a one word
!                                        integer.  Note this call does not have
!                                        an "MB" result like FMDIVI and FMMPYI.

!  FMASIN(MA,MB)        MB = ASIN(MA)

!  FMATAN(MA,MB)        MB = ATAN(MA)

!  FMATN2(MA,MB,MC)     MC = ATAN2(MA,MB)

!  FMBIG(MA)            MA = Biggest FM number less than overflow.

!  FMCHSH(MA,MB,MC)     MB = COSH(MA),  MC = SINH(MA).
!                            Faster than making two separate calls.

!  FMCOMP(MA,LREL,MB)        Logical comparison of MA and MB.
!                            LREL is a CHARACTER*2 value identifying
!                            which of the six comparisons is to be made.
!                            Example:  IF (FMCOMP(MA,'GE',MB)) ...
!                            Also can be:  IF (FMCOMP(MA,'>=',MB)) ...
!                            CHARACTER*1 is ok:  IF (FMCOMP(MA,'>',MB)) ...

!  FMCONS                    Set several saved constants that depend on MBASE,
!                            the base being used.  FMCONS should be called
!                            immediately after changing MBASE.

!  FMCOS(MA,MB)         MB = COS(MA)

!  FMCOSH(MA,MB)        MB = COSH(MA)

!  FMCSSN(MA,MB,MC)     MB = COS(MA),  MC = SIN(MA).
!                            Faster than making two separate calls.

!  FMDIG(NSTACK,KST)         Find a set of precisions to use during Newton
!                            iteration for finding a simple root starting with
!                            about double precision accuracy.

!  FMDIM(MA,MB,MC)      MC = DIM(MA,MB)

!  FMDIV(MA,MB,MC)      MC = MA / MB

!  FMDIV_R1(MA,MB)      MA = MA / MB

!  FMDIV_R2(MA,MB)      MB = MA / MB

!  FMDIVI(MA,IVAL,MB)   MB = MA/IVAL   IVAL is a one word integer.

!  FMDIVI_R1(MA,IVAL)   MA = MA/IVAL

!  FMDP2M(X,MA)         MA = X    Convert from double precision to FM.

!  FMDPM(X,MA)          MA = X    Convert from double precision to FM.
!                                 Faster than FMDP2M, but MA agrees with X only
!                                 to D.P. accuracy.  See the comments in the
!                                 two routines.

!  FMEQ(MA,MB)          MB = MA   Both have precision NDIG.
!                                 This is the version to use for standard
!                                 B = A  statements.

!  FMEQU(MA,MB,NA,NB)   MB = MA   Version for changing precision.
!                                 MA has NA digits (i.e., MA was computed
!                                 using NDIG = NA), and MB will be defined
!                                 having NB digits.
!                                 MB is rounded if NB < NA
!                                 MB is zero-padded if NB > NA

!  FMEXP(MA,MB)         MB = EXP(MA)

!  FMFLAG(K)            K = KFLAG  get the value of the FM condition
!                                  flag -- stored in the internal FM
!                                  variable KFLAG in module FMVALS.

!  FMFORM(FORM,MA,STRING)    MA is converted to a character string using format
!                               FORM and returned in STRING.  FORM can represent
!                               I, F, E, or 1PE formats.  Example:
!                               CALL FMFORM('F60.40',MA,STRING)

!  FMFPRT(FORM,MA)           Print MA on unit KW using FORM format.

!  FMI2M(IVAL,MA)       MA = IVAL   Convert from one word integer to FM.

!  FMINP(LINE,MA,LA,LB) MA = LINE   Input conversion.
!                                   Convert LINE(LA) through LINE(LB)
!                                   from characters to FM.

!  FMINT(MA,MB)         MB = INT(MA)    Integer part of MA.

!  FMIPWR(MA,IVAL,MB)   MB = MA**IVAL   Raise an FM number to a one word
!                                       integer power.

!  FMLG10(MA,MB)        MB = LOG10(MA)

!  FMLN(MA,MB)          MB = LOG(MA)

!  FMLNI(IVAL,MA)       MA = LOG(IVAL)   Natural log of a one word integer.

!  FMM2DP(MA,X)         X  = MA     Convert from FM to double precision.

!  FMM2I(MA,IVAL)       IVAL = MA   Convert from FM to integer.

!  FMM2SP(MA,X)         X  = MA     Convert from FM to single precision.

!  FMMAX(MA,MB,MC)      MC = MAX(MA,MB)

!  FMMIN(MA,MB,MC)      MC = MIN(MA,MB)

!  FMMOD(MA,MB,MC)      MC = MA mod MB

!  FMMPY(MA,MB,MC)      MC = MA * MB

!  FMMPY_R1(MA,MB)      MA = MA * MB

!  FMMPY_R2(MA,MB)      MB = MA * MB

!  FMMPYI(MA,IVAL,MB)   MB = MA*IVAL    Multiply by a one word integer.

!  FMMPYI_R1(MA,IVAL)   MA = MA*IVAL

!  FMNINT(MA,MB)        MB = NINT(MA)   Nearest FM integer.

!  FMOUT(MA,LINE,LB)    LINE = MA   Convert from FM to character.
!                                   LINE is a character array of length LB.

!  FMPI(MA)             MA = pi

!  FMPRNT(MA)                Print MA on unit KW using current format.

!  FMPWR(MA,MB,MC)      MC = MA**MB

!  FM_RANDOM_NUMBER(X)  X    is returned as a double precision random number,
!                            uniform on (0,1).  High-quality, long-period
!                            generator.
!                            Note that X is double precision, unlike the similar
!                            Fortran intrinsic random number routine, which
!                            returns a single-precision result.
!                            See the comments in section 10 below and also those
!                            in the routine for more details.

!  FMREAD(KREAD,MA)     MA   is returned after reading one (possibly multi-line)
!                            FM number on unit KREAD.  This routine reads
!                            numbers written by FMWRIT.

!  FMRPWR(MA,K,J,MB)    MB = MA**(K/J)  Rational power.
!                            Faster than FMPWR for functions like the cube root.

!  FMSET(NPREC)              Set the internal FM variables so that the precision
!                            is at least NPREC base 10 digits plus three base 10
!                            guard digits.

!  FMSETVAR(STRING)          Define a new value for one of the internal FM
!                            variables in module FMVALS that controls one of the
!                            FM options.  STRING has the form  variable = value.
!                            Example:  To change the screen width for FM output:
!                                  CALL FMSETVAR(' KSWIDE = 120 ')
!                            The variables that can be changed and the options
!                            they control are listed in sections 2 through 6
!                            above.  Only one variable can be set per call.
!                            The variable name in STRING must have no embedded
!                            blanks.  The value part of STRING can be in any
!                            numerical format, except in the case of variable
!                            CMCHAR, which is character type.  To set CMCHAR to
!                            'E', don't use any quotes in STRING:
!                                  CALL FMSETVAR(' CMCHAR = E ')

!  FMSIGN(MA,MB,MC)     MC = SIGN(MA,MB)   Sign transfer.

!  FMSIN(MA,MB)         MB = SIN(MA)

!  FMSINH(MA,MB)        MB = SINH(MA)

!  FMSP2M(X,MA)         MA = X   Convert from single precision to FM.

!  FMSQR(MA,MB)         MB = MA * MA   Faster than FMMPY.

!  FMSQR_R1(MA)         MA = MA * MA

!  FMSQRT(MA,MB)        MB = SQRT(MA)

!  FMSQRT_R1(MA)        MA = SQRT(MA)

!  FMST2M(STRING,MA)    MA = STRING
!                            Convert from character string to FM.
!                            STRING may be in any numerical format.
!                            Often more convenient than FMINP, which converts
!                            an array of CHARACTER*1 values.  Example:
!                                  CALL FMST2M('123.4',MA)

!  FMSUB(MA,MB,MC)      MC = MA - MB

!  FMSUB_R1(MA,MB)      MA = MA - MB

!  FMSUB_R2(MA,MB)      MB = MA - MB

!  FMTAN(MA,MB)         MB = TAN(MA)

!  FMTANH(MA,MB)        MB = TANH(MA)

!  FMULP(MA,MB)         MB = One Unit in the Last Place of MA.

!  FMVARS                    Write the current values of the internal FM
!                            variables on unit KW.

!  FMWRIT(KWRITE,MA)         Write MA on unit KWRITE.
!                            Multi-line numbers will have '&' as the last
!                            nonblank character on all but the last line.  These
!                            numbers can then be read easily using FMREAD.



!  These are the Gamma and Related Functions.

!  FMBERN(N,MA,MB)      MB = MA*B(N)  Multiply by Nth Bernoulli number

!  FMBETA(MA,MB,MC)     MC = Beta(MA,MB)

!  FMCOMB(MA,MB,MC)     MC = Combination MA choose MB  (Binomial coeff.)

!  FMEULR(MA)           MA = Euler's constant ( 0.5772156649... )

!  FMFACT(MA,MB)        MB = MA Factorial  (Gamma(MA+1))

!  FMGAM(MA,MB)         MB = Gamma(MA)

!  FMIBTA(MX,MA,MB,MC)  MC = Incomplete Beta(MX,MA,MB)

!  FMIGM1(MA,MB,MC)     MC = Incomplete Gamma(MA,MB).  Lower case Gamma(a,x)

!  FMIGM2(MA,MB,MC)     MC = Incomplete Gamma(MA,MB).  Upper case Gamma(a,x)

!  FMLNGM(MA,MB)        MB = Ln(Gamma(MA))

!  FMPGAM(N,MA,MB)      MB = Polygamma(N,MA)  (Nth derivative of Psi)

!  FMPOCH(MA,N,MB)      MB = MA*(MA+1)*(MA+2)*...*(MA+N-1)  (Pochhammer)

!  FMPSI(MA,MB)         MB = Psi(MA)      (Derivative of Ln(Gamma(MA))



!  These are the integer routines that are designed to be called by the user.
!  All are subroutines except logical function IMCOMP.  MA, MB, MC refer to IM
!  format numbers.  In each case the version of the routine to handle packed IM
!  numbers has the same name, with 'IM' replaced by 'IP'.

!  IMABS(MA,MB)         MB = ABS(MA)

!  IMADD(MA,MB,MC)      MC = MA + MB

!  IMBIG(MA)            MA = Biggest IM number less than overflow.

!  IMCOMP(MA,LREL,MB)        Logical comparison of MA and MB.
!                            LREL is a CHARACTER*2 value identifying which of
!                            the six comparisons is to be made.
!                            Example:  IF (IMCOMP(MA,'GE',MB)) ...
!                            Also can be:  IF (IMCOMP(MA,'>=',MB))
!                            CHARACTER*1 is ok:  IF (IMCOMP(MA,'>',MB)) ...

!  IMDIM(MA,MB,MC)      MC = DIM(MA,MB)

!  IMDIV(MA,MB,MC)      MC = int(MA/MB)
!                            Use IMDIVR if the remainder is also needed.

!  IMDIVI(MA,IVAL,MB)   MB = int(MA/IVAL)
!                            IVAL is a one word integer.
!                            Use IMDVIR to get the remainder also.

!  IMDIVR(MA,MB,MC,MD)  MC = int(MA/MB),   MD = MA mod MB
!                            When both the quotient and remainder are needed,
!                            this routine is twice as fast as calling both
!                            IMDIV and IMMOD.

!  IMDVIR(MA,IVAL,MB,IREM)   MB = int(MA/IVAL),   IREM = MA mod IVAL
!                            IVAL and IREM are one word integers.

!  IMEQ(MA,MB)          MB = MA

!  IMFM2I(MAFM,MB)      MB = MAFM  Convert from real (FM) format to
!                                  integer (IM) format.

!  IMFORM(FORM,MA,STRING)    MA is converted to a character string using format
!                               FORM and returned in STRING.  FORM can represent
!                               I, F, E, or 1PE formats.  Example:
!                               CALL IMFORM('I70',MA,STRING)

!  IMFPRT(FORM,MA)           Print MA on unit KW using FORM format.

!  IMGCD(MA,MB,MC)      MC = greatest common divisor of MA and MB.

!  IMI2FM(MA,MBFM)      MBFM = MA  Convert from integer (IM) format to
!                                  real (FM) format.

!  IMI2M(IVAL,MA)       MA = IVAL   Convert from one word integer to IM.

!  IMINP(LINE,MA,LA,LB) MA = LINE   Input conversion.
!                                   Convert LINE(LA) through LINE(LB)
!                                   from characters to IM.

!  IMM2DP(MA,X)         X  = MA     Convert from IM to double precision.

!  IMM2I(MA,IVAL)       IVAL = MA   Convert from IM to one word integer.

!  IMMAX(MA,MB,MC)      MC = MAX(MA,MB)

!  IMMIN(MA,MB,MC)      MC = MIN(MA,MB)

!  IMMOD(MA,MB,MC)      MC = MA mod MB

!  IMMPY(MA,MB,MC)      MC = MA*MB

!  IMMPYI(MA,IVAL,MB)   MB = MA*IVAL    Multiply by a one word integer.

!  IMMPYM(MA,MB,MC,MD)  MD = MA*MB mod MC
!                            Slightly faster than calling IMMPY and IMMOD
!                            separately, and it works for cases where IMMPY
!                            would return OVERFLOW.

!  IMOUT(MA,LINE,LB)    LINE = MA   Convert from IM to character.
!                                   LINE is a character array of length LB.

!  IMPMOD(MA,MB,MC,MD)  MD = MA**MB mod MC

!  IMPRNT(MA)                Print MA on unit KW.

!  IMPWR(MA,MB,MC)      MC = MA**MB

!  IMREAD(KREAD,MA)     MA   is returned after reading one (possibly multi-line)
!                            IM number on unit KREAD.
!                            This routine reads numbers written by IMWRIT.

!  IMSIGN(MA,MB,MC)     MC = SIGN(MA,MB)   Sign transfer.

!  IMSQR(MA,MB)         MB = MA*MA   Faster than IMMPY.

!  IMST2M(STRING,MA)    MA = STRING
!                            Convert from character string to IM.
!                            Often more convenient than IMINP, which converts an
!                            array of CHARACTER*1 values.  Example:
!                                 CALL IMST2M('12345678901',MA)

!  IMSUB(MA,MB,MC)      MC = MA - MB

!  IMWRIT(KWRITE,MA)         Write MA on unit KWRITE.
!                            Multi-line numbers will have '&' as the last
!                            nonblank character on all but the last line.
!                            These numbers can then be read easily using IMREAD.



!  These are the complex routines that are designed to be called by the user.
!  All are subroutines, and in each case the version of the routine to handle
!  packed ZM numbers has the same name, with 'ZM' replaced by 'ZP'.

!  MA, MB, MC refer to ZM format complex numbers.
!  MAFM, MBFM, MCFM refer to FM format real numbers.
!  INTEG is a Fortran INTEGER variable.
!  ZVAL is a Fortran COMPLEX variable.

!  ZMABS(MA,MBFM)       MBFM = ABS(MA)    Result is real.

!  ZMACOS(MA,MB)        MB = ACOS(MA)

!  ZMADD(MA,MB,MC)      MC = MA + MB

!  ZMADDI(MA,INTEG)     MA = MA + INTEG  Increment an ZM number by a one word
!                                        integer.  Note this call does not have
!                                        an "MB" result like ZMDIVI and ZMMPYI.

!  ZMARG(MA,MBFM)       MBFM = Argument(MA)    Result is real.

!  ZMASIN(MA,MB)        MB = ASIN(MA)

!  ZMATAN(MA,MB)        MB = ATAN(MA)

!  ZMCHSH(MA,MB,MC)     MB = COSH(MA),  MC = SINH(MA).
!                            Faster than 2 calls.

!  ZMCMPX(MAFM,MBFM,MC) MC = CMPLX(MAFM,MBFM)

!  ZMCONJ(MA,MB)        MB = CONJG(MA)

!  ZMCOS(MA,MB)         MB = COS(MA)

!  ZMCOSH(MA,MB)        MB = COSH(MA)

!  ZMCSSN(MA,MB,MC)     MB = COS(MA),  MC = SIN(MA).
!                            Faster than 2 calls.

!  ZMDIV(MA,MB,MC)      MC = MA / MB

!  ZMDIVI(MA,INTEG,MB)  MB = MA / INTEG

!  ZMEQ(MA,MB)          MB = MA

!  ZMEQU(MA,MB,NDA,NDB) MB = MA    Version for changing precision.
!                                  (NDA and NDB are as in FMEQU)

!  ZMEXP(MA,MB)         MB = EXP(MA)

!  ZMFORM(FORM1,FORM2,MA,STRING)   STRING = MA
!                       MA is converted to a character string using format
!                       FORM1 for the real part and FORM2 for the imaginary
!                       part.  The  result is returned in STRING.  FORM1 and
!                       FORM2 can represent I, F, E, or 1PE formats.  Example:
!                             CALL ZMFORM('F20.10','F15.10',MA,STRING)
!                       A 1PE in the first format does not carry over to the
!                       other format descriptor, as it would in an ordinary
!                       FORMAT statement.

!  ZMFPRT(FORM1,FORM2,MA)    Print MA on unit KW using formats FORM1 and FORM2.

!  ZMI2M(INTEG,MA)           MA = CMPLX(INTEG,0)

!  ZM2I2M(INTEG1,INTEG2,MA)  MA = CMPLX(INTEG1,INTEG2)

!  ZMIMAG(MA,MBFM)           MBFM = IMAG(MA)    Imaginary part.

!  ZMINP(LINE,MA,LA,LB)      MA = LINE   Input conversion.
!                                 Convert LINE(LA) through LINE(LB) from
!                                 characters to ZM.  LINE is a character array
!                                 of length at least LB.

!  ZMINT(MA,MB)         MB = INT(MA)        Integer part of both Real
!                                           and Imaginary parts of MA.

!  ZMIPWR(MA,INTEG,MB)  MB = MA ** INTEG    Integer power function.

!  ZMLG10(MA,MB)        MB = LOG10(MA)

!  ZMLN(MA,MB)          MB = LOG(MA)

!  ZMM2I(MA,INTEG)      INTEG = INT(REAL(MA))

!  ZMM2Z(MA,ZVAL)       ZVAL = MA

!  ZMMPY(MA,MB,MC)      MC = MA * MB

!  ZMMPYI(MA,INTEG,MB)  MB = MA * INTEG

!  ZMNINT(MA,MB)        MB = NINT(MA)   Nearest integer of both Real
!                                       and Imaginary.

!  ZMOUT(MA,LINE,LB,LAST1,LAST2)        LINE = MA
!                       Convert from FM to character.
!                       LINE  is the returned character*1 array.
!                       LB    is the dimensioned size of LINE.
!                       LAST1 is returned as the position in LINE of
!                             the last character of REAL(MA).
!                       LAST2 is returned as the position in LINE
!                             of the last character of AIMAG(MA).

!  ZMPRNT(MA)           Print MA on unit KW using current format.

!  ZMPWR(MA,MB,MC)      MC = MA ** MB

!  ZMREAD(KREAD,MA)     MA   is returned after reading one (possibly multi-line)
!                            ZM number on unit KREAD.
!                            This routine reads numbers written by ZMWRIT.

!  ZMREAL(MA,MBFM)      MBFM = REAL(MA)    Real part.

!  ZMRPWR(MA,IVAL,JVAL,MB)     MB = MA ** (IVAL/JVAL)

!  ZMSET(NPREC)         Set precision to the equivalent of a few more than NPREC
!                       base 10 digits.  This is now the same as FMSET, but is
!                       retained for compatibility with earlier versions of the
!                       package.

!  ZMSIN(MA,MB)         MB = SIN(MA)

!  ZMSINH(MA,MB)        MB = SINH(MA)

!  ZMSQR(MA,MB)         MB = MA*MA    Faster than ZMMPY.

!  ZMSQRT(MA,MB)        MB = SQRT(MA)

!  ZMST2M(STRING,MA)    MA = STRING
!                            Convert from character string to ZM.
!                            Often more convenient than ZMINP, which
!                            converts an array of CHARACTER*1 values.
!                            Example: CALL ZMST2M('123.4+5.67i',MA).

!  ZMSUB(MA,MB,MC)      MC = MA - MB

!  ZMTAN(MA,MB)         MB = TAN(MA)

!  ZMTANH(MA,MB)        MB = TANH(MA)

!  ZMWRIT(KWRITE,MA)    Write MA on unit KWRITE.  Multi-line numbers are
!                       formatted for automatic reading with ZMREAD.

!  ZMZ2M(ZVAL,MA)       MA = ZVAL


!  10. NEW FOR VERSION 1.2

!  Version 1.2 is written in Fortran-90 free source format.

!  The routines for the Gamma function and related mathematical special
!  functions are new in version 1.2.

!  Several new derived-type function interfaces are included in module FMZM in
!  file FMZM90.f90, such as integer multiple precision operations GCD, modular
!  multiplication, and modular powers.  There are also formatting functions and
!  function interfaces for the Gamma and related special functions.

!  Two new rounding modes have been added, round toward -infinity and round
!  toward +infinity.  See the description of KROUND above.
!  An option has been added to force more guard digits to be used, so that basic
!  arithmetic operations will always round perfectly.  See the description of
!  KRPERF above.
!  These options are included for applications that use FM to check IEEE
!  hardware arithmetic.  They are not normally useful for most multiple
!  precision calculations.

!  The random number routine FM_RANDOM_NUMBER uses 49-digit prime numbers in a
!  shuffled multiplicative congruential generator.  Historically, some popular
!  random number routines tried so hard for maximum speed that they were later
!  found to fail some tests for randomness.  FM_RANDOM_NUMBER tries to return
!  high-quality random values.  It is much slower than other generators, but can
!  return about 60,000 numbers per second on a 400 MHz single-processor machine.
!  This is usually fast enough to be used as a check for suspicious monte carlo
!  results from other generators.
!  For more details, see the comments in the routine.

!  The arrays for multiple precision numbers were dimensioned starting at 0 in
!  version 1.1, and now begin at -1.  Array(-1) now holds the sign of the number
!  instead of combining the sign with Array(2) as before.  The reason for moving
!  the sign bit is that many of the original routines, written before Fortran-90
!  existed, simplified the logic by temporarily making input arguments positive,
!  working with positive values, then restoring the signs to the input arguments
!  upon return.  This became illegal under Fortran-90 when used with the derived
!  type interface, which demands the inputs to functions for arithmetic operator
!  overloading be declared with INTENT(IN).

!  The common blocks of earlier versions have been replaced by module FMVALS.
!  This makes it easier to hide the FM internal variable names from the calling
!  program, and these variables can be initialized in the module so the
!  initializing call to FMSET is no longer mandatory.  Several new routines are
!  provided to set or return the values for some of these variables.  See the
!  descriptions for FMSETVAR, FMFLAG, and FMVARS above.

!  Version 1.0 used integer arrays and integer arithmetic internally to perform
!  the multiple precision operations.  Later versions use double precision
!  arithmetic and arrays internally.  This is usually faster at higher
!  precisions, and on many machines it is also faster at lower precisions.
!  Version 1.2 is written so that the arithmetic used can easily be changed from
!  double precision to integer, or any other available arithmetic type.  This
!  permits the user to make the best use of a given machine's arithmetic
!  hardware.  See the EFFICIENCY discussion below.


!  11. EFFICIENCY

!  When the derived type interface is used to access the FM routines, there may
!  be a loss of speed if the arrays used to define the multiple precision data
!  types are larger than necessary.  See comment (b) in the section above on
!  array dimensions.

!  To take advantage of hardware architecture on different machines, the package
!  has been designed so that the arithmetic used to perform the multiple
!  precision operations can easily be changed.  All variables that must be
!  changed to get a different arithmetic have names beginning with 'M' and are
!  declared using REAL (KIND(1.0D0)) ...

!  For example, to change the package to use integer arithmetic internally, make
!  these two changes everywhere in the FM.f90 file.
!  Change  'REAL (KIND(1.0D0))'  to  'INTEGER'.
!  Change  'AINT ('  to  'INT('.  Note the blank between AINT and (.
!  On some systems, changing  'AINT ('  to  '('  may give better speed.

!  In most places in FM, an AINT function is not supposed to be changed.  These
!  are written 'AINT(', with no embedded blank, so they will not be changed by
!  the global change above.

!  The first of these changes must also be made throughout the files FMZM90.f90
!  and FMSAVE.f90.
!  Change  'REAL (KIND(1.0D0))'  to  'INTEGER'.

!  Many of the variables in FMSAVE.f90 are initialized when they are declared,
!  so the initialization values should be changed to integer values.
!  Find the lines beginning '! Integer initialization' in file FMSAVE.f90 and
!  change the values.  The values needed for 32-bit integer arithmetic are next
!  to the double precision values, but commented out.  In every case, the line
!  before the '! Integer initialization' should have '!' inserted in column 1
!  and the line after should have the '!' removed from column 1.  If a different
!  wordsize is used, the first call to FMSET will check the values defined in
!  file FMSAVE.f90 and write messages (on unit KW) if any need to be changed.

!  When changing to a different type of arithmetic, any FM arrays in the user's
!  program must be changed to agree.  If derived types are used instead of
!  direct calls, no changes should be needed in the calling program.

!  For example, in the test program TestFM.f90, change all
!  'REAL (KIND(1.0D0))'  to  'INTEGER', as with the other files.

!  This version of FM restricts the base used to be also representable in
!  integer variables, so using precision above double usually does not save much
!  time unless integers can also be declared at a higher precision.  Using IEEE
!  Extended would allow a base of around 10**9 to be chosen, but the delayed
!  digit-normalization method used for multiplication and division means that a
!  slightly smaller base like 10**8 would probably run faster.  This would
!  usually not be much faster than using the usual base 10**7 with double
!  precision.

!  The value of NBITS defined as a parameter in FMVALS refers to the number of
!  bits used to represent integers in an M-variable word.  Typical values for
!  NBITS are:  24 for IEEE single precision, 32 for integer, 53 for IEEE double
!  precision.  NBITS controls only array size, so setting it too high is ok, but
!  then the program will use slightly more memory than necessary.

!  For cases where special compiler directives or minor re-writing of the code
!  may improve speed, several of the most important loops in FM are identified
!  by comments containing the string '(Inner Loop)'.

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


      SUBROUTINE FMSET(NPREC)

!  Initialize the global FM variables that must be set before calling
!  other FM routines.  These variables are initialized to fairly standard
!  values in the FMSAVE.f90 file (MODULE FMVALS), so calling FMSET at the
!  beginning of a program is now optional.  FMSET is a convenient way to set
!  or change the precision being used, and it also checks to see that the
!  generic values chosen for several machine-dependent variables are valid.

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

!  MBASE (base for FM arithmetic) is set to a large power of ten.
!  JFORM1 and JFORM2 (default output format controls) are set to 1PE format
!  displaying NPREC significant digits.

!  Several FM options were set here in previous versions of the package,
!  and are now initialized to their default values in module FMVALS.
!  Here are the initial settings:

!  The trace option is set off.
!  The mode for angles in trig functions is set to radians.
!  The rounding mode is set to symmetric rounding.
!  Warning error message level is set to 1.
!  Cancellation error monitor is set off.
!  Screen width for output is set to 80 columns.
!  The exponent character for FM output is set to 'M'.
!  Debug error checking is set off.

      USE FMVALS
      IMPLICIT NONE

      INTEGER NPREC

      REAL (KIND(1.0D0)) :: MAXINT_CHK,MXEXP2_CHK,MEXPOV_CHK,MEXPUN_CHK,  &
                            MUNKNO_CHK
      DOUBLE PRECISION DPEPS_CHK,DPMAX_CHK,SPMAX_CHK,TEMP
      INTEGER INTMAX_CHK,K,NPSAVE

      IF (NBITS < DIGITS(MAXINT)) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' NBITS was set to ',NBITS,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be at least ',DIGITS(MAXINT)
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' '
          WRITE (KW,*) ' NBITS is a parameter that controls array size, so its'
          WRITE (KW,*) ' value cannot be changed for this run, and this might'
          WRITE (KW,*) ' cause some FM operations to get incorrect results.'
          WRITE (KW,*) ' '
      ENDIF

!             MAXINT should be set to a very large integer, possibly
!                    the largest representable integer for the current
!                    machine.  For most 32-bit machines, MAXINT is set
!                    to  2**53 - 1 = 9.007D+15  when double precision
!                    arithmetic is used for M-variables.  Using integer
!                    M-variables usually gives MAXINT = 2**31 - 1 =
!                    2147483647.

!                    Setting MAXINT to a smaller number is ok, but this
!                    unnecessarily restricts the permissible range of
!                    MBASE and MXEXP.

      MAXINT_CHK = RADIX(MAXINT_CHK)
      MAXINT_CHK = ((MAXINT_CHK**(DIGITS(MAXINT_CHK)-1)-1)*MAXINT_CHK - 1) +  &
                     MAXINT_CHK
      IF (MAXINT > MAXINT_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MAXINT was set to ',MAXINT,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',MAXINT_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MAXINT has been changed to ',MAXINT_CHK
          WRITE (KW,*) ' '
          MAXINT = MAXINT_CHK
      ELSE IF (MAXINT < MAXINT_CHK/2) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MAXINT was set to ',MAXINT,' in file FMSAVE.f90'
          WRITE (KW,*) ' For better performance set it to ',MAXINT_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MAXINT has been changed to ',MAXINT_CHK
          WRITE (KW,*) ' '
          MAXINT = MAXINT_CHK
      ENDIF

!             INTMAX is a large value close to the overflow threshold
!                    for integer variables.  It is usually 2**31 - 1
!                    for machines with 32-bit integer arithmetic.

!                    The following code sets INTMAX_CHK to the
!                    largest representable integer.
!                    Then INTMAX is checked against this value.

      INTMAX_CHK = HUGE(1)
      IF (INTMAX > INTMAX_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' INTMAX was set to ',INTMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',INTMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, INTMAX has been changed to ',INTMAX_CHK
          WRITE (KW,*) ' '
          INTMAX = INTMAX_CHK
      ELSE IF (INTMAX < INTMAX_CHK/2) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' INTMAX was set to ',INTMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For better performance set it to ',INTMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, INTMAX has been changed to ',INTMAX_CHK
          WRITE (KW,*) ' '
          INTMAX = INTMAX_CHK
      ENDIF

!             DPMAX should be set to a value near the machine's double
!                   precision overflow threshold, so that DPMAX and
!                   1.0D0/DPMAX are both representable in double
!                   precision.

      DPMAX_CHK = HUGE(1.0D0)/5
      IF (DPMAX > DPMAX_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' DPMAX was set to ',DPMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',DPMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, DPMAX has been changed to ',DPMAX_CHK
          WRITE (KW,*) ' '
          DPMAX = DPMAX_CHK
      ELSE IF (DPMAX < DPMAX_CHK/1.0D2) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' DPMAX was set to ',DPMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For better performance set it to ',DPMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, DPMAX has been changed to ',DPMAX_CHK
          WRITE (KW,*) ' '
          DPMAX = DPMAX_CHK
      ENDIF

!             SPMAX should be set to a value near the machine's single
!                   precision overflow threshold, so that 1.01*SPMAX
!                   and 1.0/SPMAX are both representable in single
!                   precision.

      SPMAX_CHK = HUGE(1.0)/5
      IF (SPMAX > SPMAX_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' SPMAX was set to ',SPMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',SPMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, SPMAX has been changed to ',SPMAX_CHK
          WRITE (KW,*) ' '
          SPMAX = SPMAX_CHK
      ELSE IF (SPMAX < SPMAX_CHK/1.0D2) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' SPMAX was set to ',SPMAX,' in file FMSAVE.f90'
          WRITE (KW,*) ' For better performance set it to ',SPMAX_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, SPMAX has been changed to ',SPMAX_CHK
          WRITE (KW,*) ' '
          SPMAX = SPMAX_CHK
      ENDIF

!             MXBASE is the maximum value for MBASE.

      TEMP = MAXINT
      TEMP = INT(MIN(DBLE(INTMAX),SQRT(TEMP)))
      IF (MXBASE > TEMP) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MXBASE was set to ',MXBASE,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',TEMP
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MXBASE has been changed to ',TEMP
          WRITE (KW,*) ' '
          MXBASE = TEMP
      ELSE IF (MXBASE < TEMP/2) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MXBASE was set to ',MXBASE,' in file FMSAVE.f90'
          WRITE (KW,*) ' For better performance set it to ',TEMP
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MXBASE has been changed to ',TEMP
          WRITE (KW,*) ' '
          MXBASE = TEMP
      ENDIF

!             MBASE is the currently used base for arithmetic.

      K = INT(LOG10(DBLE(MXBASE)/4))
      MBASE = 10**K

!             NDIG is the number of digits currently being carried.

      NPSAVE = NPREC
      NDIG = 2 + (NPREC+2)/K
      IF (NDIG < 2 .OR. NDIG > NDIGMX) THEN
          NDIG = MAX(2,NDIG)
          NDIG = MIN(NDIGMX,NDIG)
          WRITE (KW,  &
                 "(//' Precision out of range when calling FMSET.',"  //  &
                 "'  NPREC =',I20/' The nearest valid NDIG will be'," //  &
                 "' used instead:   NDIG =',I6//)"                        &
                ) NPREC,NDIG
          NPSAVE = 0
      ENDIF

!             NCALL is the call stack pointer.

      NCALL = 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(INTMAX) is greater
!             than MXBASE**(MXEXP2+1).

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

      MXEXP = INT((DBLE(INTMAX))/(2.0D0*LOG(DBLE(MXBASE))) - 1.0D0)
      MXEXP2_CHK = INT(2*MXEXP + MXEXP/100)
      IF (MXEXP2 > MXEXP2_CHK*1.01) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MXEXP2 was set to ',MXEXP2,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',MXEXP2_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MXEXP2 has been changed to ',MXEXP2_CHK
          WRITE (KW,*) ' '
          MXEXP2 = MXEXP2_CHK
      ELSE IF (MXEXP2 < MXEXP2_CHK*0.99) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MXEXP2 was set to ',MXEXP2,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no less than ',MXEXP2_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MXEXP2 has been changed to ',MXEXP2_CHK
          WRITE (KW,*) ' '
          MXEXP2 = MXEXP2_CHK
      ENDIF

!             KACCSW is a switch used to enable cancellation error
!                    monitoring.  Routines where cancellation is
!                    not a problem run faster by skipping the
!                    cancellation monitor calculations.
!                    KACCSW = 0 means no error monitoring,
!                           = 1 means error monitoring is done.

      KACCSW = 0

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

      MEXPUN_CHK = -AINT(MXEXP2*1.01D0)
      IF (MEXPUN < MEXPUN_CHK*1.01) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MEXPUN was set to ',MEXPUN,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no less than ',MEXPUN_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MEXPUN has been changed to ',MEXPUN_CHK
          WRITE (KW,*) ' '
          MEXPUN = MEXPUN_CHK
      ELSE IF (MEXPUN > MEXPUN_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MEXPUN was set to ',MEXPUN,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',MEXPUN_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MEXPUN has been changed to ',MEXPUN_CHK
          WRITE (KW,*) ' '
          MEXPUN = MEXPUN_CHK
      ENDIF

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

      MEXPOV_CHK = -MEXPUN
      IF (MEXPOV /= MEXPOV_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MEXPOV was set to ',MEXPOV,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be ',MEXPOV_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MEXPOV has been changed to ',MEXPOV_CHK
          WRITE (KW,*) ' '
          MEXPOV = MEXPOV_CHK
      ENDIF

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

      MUNKNO_CHK = AINT(MEXPOV*1.01D0)
      IF (MUNKNO > MUNKNO_CHK*1.01) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MUNKNO was set to ',MUNKNO,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',MUNKNO_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MUNKNO has been changed to ',MUNKNO_CHK
          WRITE (KW,*) ' '
          MUNKNO = MUNKNO_CHK
      ELSE IF (MUNKNO < MUNKNO_CHK) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' MUNKNO was set to ',MUNKNO,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no less than ',MUNKNO_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, MUNKNO has been changed to ',MUNKNO_CHK
          WRITE (KW,*) ' '
          MUNKNO = MUNKNO_CHK
      ENDIF

!             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 that 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 that is likely to
!                    make any subsequent results obviously wrong that
!                    use it.  In either case a KFLAG = -4 condition is
!                    also returned.

      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.

      IUNKNO = -INT(MXEXP2)

!             DPEPS is the approximate machine precision.

      DPEPS_CHK = EPSILON(1.0D0)
      IF (DPEPS > DPEPS_CHK*1.01) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' DPEPS was set to ',DPEPS,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no more than ',DPEPS_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, DPEPS has been changed to ',DPEPS_CHK
          WRITE (KW,*) ' '
          DPEPS = DPEPS_CHK
      ELSE IF (DPEPS < DPEPS_CHK*0.99) THEN
          WRITE (KW,*) ' '
          WRITE (KW,*) ' In routine FMSET it appears that FM internal variable'
          WRITE (KW,*) ' DPEPS was set to ',DPEPS,' in file FMSAVE.f90'
          WRITE (KW,*) ' For this machine it should be no less than ',DPEPS_CHK
          WRITE (KW,*) ' Change the initialization in FMSAVE.f90 to this value.'
          WRITE (KW,*) ' For this run, DPEPS has been changed to ',DPEPS_CHK
          WRITE (KW,*) ' '
          DPEPS = DPEPS_CHK
      ENDIF

!             JFORM1 indicates the format used by FMOUT.

      JFORM1 = 1

!             JFORM2 indicates the number of digits used in FMOUT.

      JFORM2 = NPSAVE

!             Set JFORMZ to ' 1.23 + 4.56 i ' format.

      JFORMZ = 1

!             Set JPRNTZ to print real and imaginary parts on one
!             line whenever possible.

      JPRNTZ = 1

!             Initialize two hash tables that are used for character
!             look-up during input conversion.

      CALL FMHTBL

!             FMCONS sets several real and double precision constants.

      CALL FMCONS


      RETURN
      END SUBROUTINE FMSET

      SUBROUTINE FMABS(MA,MB)

!  MB = ABS(MA)

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)

      REAL (KIND(1.0D0)) :: MD2B
      INTEGER KWRNSV

      NCALL = NCALL + 1
      NAMEST(NCALL) = 'FMABS '
      IF (NTRACE /= 0) CALL FMNTR(2,MA,MA,1,1)

      KFLAG = 0
      KWRNSV = KWARN
      KWARN = 0
      CALL FMEQ(MA,MB)
      MB(-1) = 1
      KWARN = KWRNSV

      IF (KACCSW == 1) THEN
          MD2B = NINT((NDIG-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)
          MB(0) = MIN(MB(0),MD2B)
      ENDIF
      IF (NTRACE /= 0) CALL FMNTR(1,MB,MB,1,1)
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMABS

      SUBROUTINE FMACOS(MA,MB)

!  MB = ACOS(MA)

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      REAL (KIND(1.0D0)) :: MACCA,MACMAX,MAS,MXSAVE
      INTEGER J,K,KASAVE,KOVUN,KRESLT,NDSAVE

      IF (MBLOGS /= MBASE) CALL FMCONS
      IF (ABS(MA(1)) > MEXPAB .OR. MA(1) > 0 .OR. MA(2) == 0) THEN
          CALL FMENTR('FMACOS',MA,MA,1,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,  &
                      KOVUN)
          IF (KRESLT /= 0) RETURN
      ELSE
          NCALL = NCALL + 1
          NAMEST(NCALL) = 'FMACOS'
          IF (NTRACE /= 0) CALL FMNTR(2,MA,MA,1,1)
          KOVUN = 0
          IF (MA(1) == MEXPOV .OR. MA(1) == MEXPUN) KOVUN = 1
          NDSAVE = NDIG
          IF (NCALL == 1) THEN
              K = MAX(NGRD52-1,2)
              NDIG = MAX(NDIG+K,2)
              IF (NDIG > NDG2MX) THEN
                  KFLAG = -9
                  CALL FMWARN
                  NDIG = NDSAVE
                  KRESLT = 12
                  CALL FMRSLT(MA,MA,MB,KRESLT)
                  IF (NTRACE /= 0) CALL FMNTR(1,MB,MB,1,1)
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
          KASAVE = KACCSW
          KACCSW = 0
          MXSAVE = MXEXP
          MXEXP = MXEXP2
      ENDIF

      MAS = MA(-1)
      MACCA = MA(0)
      CALL FMEQ2(MA,MB,NDSAVE,NDIG)
      MB(0) = NINT(NDIG*ALOGM2)

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

      MB(-1) = 1
      CALL FMI2M(1,M05)
      CALL FMSUB(M05,MB,M03)
      CALL FMADD(M05,MB,M04)
      CALL FMMPY_R2(M03,M04)
      CALL FMSQRT_R1(M04)
      CALL FMDIV_R2(M04,MB)

      CALL FMATAN(MB,M13)
      CALL FMEQ(M13,MB)

      IF (MAS < 0) THEN
          IF (KRAD == 1) THEN
              CALL FMPI(M05)
          ELSE
              CALL FMI2M(180,M05)
          ENDIF
          CALL FMSUB_R2(M05,MB)
      ENDIF

!             Round the result and return.

      MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)
      MB(0) = MIN(MB(0),MACCA,MACMAX)
      DO J = -1, NDIG+1
         M01(J) = MB(J)
      ENDDO
      CALL FMEXIT(M01,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
      RETURN
      END SUBROUTINE FMACOS

      SUBROUTINE FMADD(MA,MB,MC)

!  MC = MA + MB

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

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK),MC(-1:LUNPCK)

      NCALL = NCALL + 1
      IF (NTRACE /= 0) THEN
          NAMEST(NCALL) = 'FMADD '
          CALL FMNTR(2,MA,MB,2,1)

          CALL FMADD2(MA,MB,MC)

          CALL FMNTR(1,MC,MC,1,1)
      ELSE
          CALL FMADD2(MA,MB,MC)
      ENDIF
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMADD

      SUBROUTINE FMADD2(MA,MB,MC)

!  Internal addition routine.  MC = MA + MB

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK),MC(-1:LUNPCK)

      REAL (KIND(1.0D0)) :: MA0,MA1,MA2,MAS,MB0,MB1,MB2,MB2RD,MBS
      INTEGER J,JCOMP,JRSSAV,JSIGN,KRESLT,N1,NGUARD,NMWA
      REAL B2RDA,B2RDB

      IF (MBLOGS /= MBASE) CALL FMCONS
      JRSSAV = JRSIGN
      MA2 = MA(2)
      MB2 = MB(2)
      IF (ABS(MA(1)) > MEXPAB .OR. ABS(MB(1)) > MEXPAB .OR.  &
          KDEBUG == 1) THEN
          IF (KSUB == 1) THEN
              CALL FMARGS('FMSUB ',2,MA,MB,KRESLT)
          ELSE
              CALL FMARGS('FMADD ',2,MA,MB,KRESLT)
          ENDIF
          IF (KRESLT /= 0) THEN
              IF ((KRESLT /= 1 .AND. KRESLT /= 2) .OR. MA(2) == 0 .OR.  &
                  MB(2) == 0) THEN
                  NCALL = NCALL + 1
                  IF (KSUB == 1) THEN
                      NAMEST(NCALL) = 'FMSUB '
                  ELSE
                      NAMEST(NCALL) = 'FMADD '
                  ENDIF
                  CALL FMRSLT(MA,MB,MC,KRESLT)
                  JRSIGN = JRSSAV
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
      ELSE
          IF (MA(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              CALL FMEQ(MB,MC)
              MC(0) = MA0
              KFLAG = 1
              IF (KSUB == 1) THEN
                  IF (MC(1) /= MUNKNO .AND. MC(2) /= 0) MC(-1) = -MC(-1)
                  KFLAG = 0
              ENDIF
              JRSIGN = JRSSAV
              RETURN
          ENDIF
          IF (MB(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              CALL FMEQ(MA,MC)
              MC(0) = MA0
              KFLAG = 1
              JRSIGN = JRSSAV
              RETURN
          ENDIF
      ENDIF

      MA0 = MA(0)
      IF (KACCSW == 1) THEN
          MB0 = MB(0)
          MA1 = MA(1)
          MB1 = MB(1)
      ENDIF
      KFLAG = 0
      N1 = NDIG + 1

!             NGUARD is the number of guard digits used.

      IF (NCALL > 1) THEN
          NGUARD = NGRD21
          IF (NGUARD > NDIG) NGUARD = NDIG
      ELSE
          NGUARD = NGRD52
          IF (NGUARD > NDIG) NGUARD = NDIG
          IF (MBASE < 1000 .OR. KROUND /= 1 .OR. KRPERF == 1) THEN
              NGUARD = NDIG + 10
          ENDIF
      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
      MAS = MA(-1)
      MBS = MB(-1)
      IF (KSUB == 1) MBS = -MBS

!             See which one is larger in absolute value.

      JCOMP = 2
      IF (MA(1) > MB(1)) THEN
          JCOMP = 1
      ELSE IF (MB(1) > MA(1)) THEN
          JCOMP = 3
      ELSE
          DO J = 2, N1
             IF (MA(J) > MB(J)) THEN
                 JCOMP = 1
                 EXIT
             ENDIF
             IF (MB(J) > MA(J)) THEN
                 JCOMP = 3
                 EXIT
             ENDIF
          ENDDO
      ENDIF

      IF (JCOMP < 3) THEN
          IF (MAS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MA,MB,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MA,MB,NGUARD,NMWA)
          ENDIF
      ELSE
          IF (MBS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MB,MA,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MB,MA,NGUARD,NMWA)
          ENDIF
      ENDIF

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

      CALL FMMOVE(MWA,MC)
      MC(-1) = 1
      IF (JSIGN < 0 .AND. MC(2) /= 0) MC(-1) = -1

      IF (KFLAG < 0) THEN
          IF (KSUB == 1) THEN
              NAMEST(NCALL) = 'FMSUB '
          ELSE
              NAMEST(NCALL) = 'FMADD '
          ENDIF
          CALL FMWARN
      ENDIF

      IF (KACCSW == 1) THEN
          B2RDA = LOG(REAL(ABS(MC(2))+1)/REAL(ABS(MA2)+1))/0.69315 +  &
                  REAL(MC(1)-MA1)*ALOGM2 + REAL(MA0)
          B2RDB = LOG(REAL(ABS(MC(2))+1)/REAL(ABS(MB2)+1))/0.69315 +  &
                  REAL(MC(1)-MB1)*ALOGM2 + REAL(MB0)
          MB2RD = NINT(MAX(0.0,MIN(B2RDA,B2RDB,  &
                  (NDIG-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315)))
          IF (MC(2) == 0) THEN
              MC(0) = 0
          ELSE
              MC(0) = MIN(MAX(MA0,MB0),MB2RD)
          ENDIF
      ELSE
          MC(0) = MA0
      ENDIF

      JRSIGN = JRSSAV
      RETURN
      END SUBROUTINE FMADD2

      SUBROUTINE FMADD_R1(MA,MB)

!  MA = MA + MB

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

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)

      NCALL = NCALL + 1
      IF (NTRACE /= 0) THEN
          NAMEST(NCALL) = 'FMADD '
          CALL FMNTR(2,MA,MB,2,1)

          CALL FMADD2_R1(MA,MB)

          CALL FMNTR(1,MA,MA,1,1)
      ELSE
          CALL FMADD2_R1(MA,MB)
      ENDIF
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMADD_R1

      SUBROUTINE FMADD2_R1(MA,MB)

!  Internal addition routine.  MA = MA + MB

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)

      REAL (KIND(1.0D0)) :: MA0,MA1,MA2,MAS,MB0,MB1,MB2,MB2RD,MBS
      INTEGER J,JCOMP,JRSSAV,JSIGN,KRESLT,N1,NGUARD,NMWA
      REAL B2RDA,B2RDB

      IF (MBLOGS /= MBASE) CALL FMCONS
      JRSSAV = JRSIGN
      MA2 = MA(2)
      MB2 = MB(2)
      IF (ABS(MA(1)) > MEXPAB .OR. ABS(MB(1)) > MEXPAB .OR.  &
          KDEBUG == 1) THEN
          IF (KSUB == 1) THEN
              CALL FMARGS('FMSUB ',2,MA,MB,KRESLT)
          ELSE
              CALL FMARGS('FMADD ',2,MA,MB,KRESLT)
          ENDIF
          IF (KRESLT /= 0) THEN
              IF ((KRESLT /= 1 .AND. KRESLT /= 2) .OR. MA(2) == 0 .OR.  &
                  MB(2) == 0) THEN
                  NCALL = NCALL + 1
                  IF (KSUB == 1) THEN
                      NAMEST(NCALL) = 'FMSUB '
                  ELSE
                      NAMEST(NCALL) = 'FMADD '
                  ENDIF
                  CALL FMRSLT(MA,MB,M07,KRESLT)
                  CALL FMEQ(M07,MA)
                  JRSIGN = JRSSAV
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
      ELSE
          IF (MA(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              CALL FMEQ(MB,MA)
              MA(0) = MA0
              KFLAG = 1
              IF (KSUB == 1) THEN
                  IF (MA(1) /= MUNKNO .AND. MA(2) /= 0) MA(-1) = -MA(-1)
                  KFLAG = 0
              ENDIF
              JRSIGN = JRSSAV
              RETURN
          ENDIF
          IF (MB(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              MA(0) = MA0
              KFLAG = 1
              JRSIGN = JRSSAV
              RETURN
          ENDIF
      ENDIF

      MA0 = MA(0)
      IF (KACCSW == 1) THEN
          MB0 = MB(0)
          MA1 = MA(1)
          MB1 = MB(1)
      ENDIF
      KFLAG = 0
      N1 = NDIG + 1

!             NGUARD is the number of guard digits used.

      IF (NCALL > 1) THEN
          NGUARD = NGRD21
          IF (NGUARD > NDIG) NGUARD = NDIG
      ELSE
          NGUARD = NGRD52
          IF (NGUARD > NDIG) NGUARD = NDIG
          IF (MBASE < 1000 .OR. KROUND /= 1 .OR. KRPERF == 1) THEN
              NGUARD = NDIG + 10
          ENDIF
      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
      MAS = MA(-1)
      MBS = MB(-1)
      IF (KSUB == 1) MBS = -MBS

!             See which one is larger in absolute value.

      JCOMP = 2
      IF (MA(1) > MB(1)) THEN
          JCOMP = 1
      ELSE IF (MB(1) > MA(1)) THEN
          JCOMP = 3
      ELSE
          DO J = 2, N1
             IF (MA(J) > MB(J)) THEN
                 JCOMP = 1
                 EXIT
             ENDIF
             IF (MB(J) > MA(J)) THEN
                 JCOMP = 3
                 EXIT
             ENDIF
          ENDDO
      ENDIF

      IF (JCOMP < 3) THEN
          IF (MAS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MA,MB,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MA,MB,NGUARD,NMWA)
          ENDIF
      ELSE
          IF (MBS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MB,MA,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MB,MA,NGUARD,NMWA)
          ENDIF
      ENDIF

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

      CALL FMMOVE(MWA,MA)
      MA(-1) = 1
      IF (JSIGN < 0 .AND. MA(2) /= 0) MA(-1) = -1

      IF (KFLAG < 0) THEN
          IF (KSUB == 1) THEN
              NAMEST(NCALL) = 'FMSUB '
          ELSE
              NAMEST(NCALL) = 'FMADD '
          ENDIF
          CALL FMWARN
      ENDIF

      IF (KACCSW == 1) THEN
          B2RDA = LOG(REAL(ABS(MA(2))+1)/REAL(ABS(MA2)+1))/0.69315 +  &
                  REAL(MA(1)-MA1)*ALOGM2 + REAL(MA0)
          B2RDB = LOG(REAL(ABS(MA(2))+1)/REAL(ABS(MB2)+1))/0.69315 +  &
                  REAL(MA(1)-MB1)*ALOGM2 + REAL(MB0)
          MB2RD = NINT(MAX(0.0,MIN(B2RDA,B2RDB,  &
                  (NDIG-1)*ALOGM2 + LOG(REAL(ABS(MA(2))+1))/0.69315)))
          IF (MA(2) == 0) THEN
              MA(0) = 0
          ELSE
              MA(0) = MIN(MAX(MA0,MB0),MB2RD)
          ENDIF
      ELSE
          MA(0) = MA0
      ENDIF

      JRSIGN = JRSSAV
      RETURN
      END SUBROUTINE FMADD2_R1

      SUBROUTINE FMADD_R2(MA,MB)

!  MB = MA + MB

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

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)

      NCALL = NCALL + 1
      IF (NTRACE /= 0) THEN
          NAMEST(NCALL) = 'FMADD '
          CALL FMNTR(2,MA,MB,2,1)

          CALL FMADD2_R2(MA,MB)

          CALL FMNTR(1,MB,MB,1,1)
      ELSE
          CALL FMADD2_R2(MA,MB)
      ENDIF
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMADD_R2

      SUBROUTINE FMADD2_R2(MA,MB)

!  Internal addition routine.  MB = MA + MB

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)

      REAL (KIND(1.0D0)) :: MA0,MA1,MA2,MAS,MB0,MB1,MB2,MB2RD,MBS
      INTEGER J,JCOMP,JRSSAV,JSIGN,KRESLT,N1,NGUARD,NMWA
      REAL B2RDA,B2RDB

      IF (MBLOGS /= MBASE) CALL FMCONS
      JRSSAV = JRSIGN
      MA2 = MA(2)
      MB2 = MB(2)
      IF (ABS(MA(1)) > MEXPAB .OR. ABS(MB(1)) > MEXPAB .OR.  &
          KDEBUG == 1) THEN
          IF (KSUB == 1) THEN
              CALL FMARGS('FMSUB ',2,MA,MB,KRESLT)
          ELSE
              CALL FMARGS('FMADD ',2,MA,MB,KRESLT)
          ENDIF
          IF (KRESLT /= 0) THEN
              IF ((KRESLT /= 1 .AND. KRESLT /= 2) .OR. MA(2) == 0 .OR.  &
                  MB(2) == 0) THEN
                  NCALL = NCALL + 1
                  IF (KSUB == 1) THEN
                      NAMEST(NCALL) = 'FMSUB '
                  ELSE
                      NAMEST(NCALL) = 'FMADD '
                  ENDIF
                  CALL FMRSLT(MA,MB,M07,KRESLT)
                  CALL FMEQ(M07,MB)
                  JRSIGN = JRSSAV
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
      ELSE
          IF (MA(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              MB(0) = MA0
              KFLAG = 1
              IF (KSUB == 1) THEN
                  IF (MB(1) /= MUNKNO .AND. MB(2) /= 0) MB(-1) = -MB(-1)
                  KFLAG = 0
              ENDIF
              JRSIGN = JRSSAV
              RETURN
          ENDIF
          IF (MB(2) == 0) THEN
              MA0 = MIN(MA(0),MB(0))
              CALL FMEQ(MA,MB)
              MB(0) = MA0
              KFLAG = 1
              JRSIGN = JRSSAV
              RETURN
          ENDIF
      ENDIF

      MA0 = MA(0)
      IF (KACCSW == 1) THEN
          MB0 = MB(0)
          MA1 = MA(1)
          MB1 = MB(1)
      ENDIF
      KFLAG = 0
      N1 = NDIG + 1

!             NGUARD is the number of guard digits used.

      IF (NCALL > 1) THEN
          NGUARD = NGRD21
          IF (NGUARD > NDIG) NGUARD = NDIG
      ELSE
          NGUARD = NGRD52
          IF (NGUARD > NDIG) NGUARD = NDIG
          IF (MBASE < 1000 .OR. KROUND /= 1 .OR. KRPERF == 1) THEN
              NGUARD = NDIG + 10
          ENDIF
      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
      MAS = MA(-1)
      MBS = MB(-1)
      IF (KSUB == 1) MBS = -MBS

!             See which one is larger in absolute value.

      JCOMP = 2
      IF (MA(1) > MB(1)) THEN
          JCOMP = 1
      ELSE IF (MB(1) > MA(1)) THEN
          JCOMP = 3
      ELSE
          DO J = 2, N1
             IF (MA(J) > MB(J)) THEN
                 JCOMP = 1
                 EXIT
             ENDIF
             IF (MB(J) > MA(J)) THEN
                 JCOMP = 3
                 EXIT
             ENDIF
          ENDDO
      ENDIF

      IF (JCOMP < 3) THEN
          IF (MAS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MA,MB,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MA,MB,NGUARD,NMWA)
          ENDIF
      ELSE
          IF (MBS < 0) JSIGN = -1
          JRSIGN = JSIGN
          IF (MAS*MBS > 0) THEN
              CALL FMADDP(MB,MA,NGUARD,NMWA)
          ELSE
              CALL FMADDN(MB,MA,NGUARD,NMWA)
          ENDIF
      ENDIF

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

      CALL FMMOVE(MWA,MB)
      MB(-1) = 1
      IF (JSIGN < 0 .AND. MB(2) /= 0) MB(-1) = -1

      IF (KFLAG < 0) THEN
          IF (KSUB == 1) THEN
              NAMEST(NCALL) = 'FMSUB '
          ELSE
              NAMEST(NCALL) = 'FMADD '
          ENDIF
          CALL FMWARN
      ENDIF

      IF (KACCSW == 1) THEN
          B2RDA = LOG(REAL(ABS(MB(2))+1)/REAL(ABS(MA2)+1))/0.69315 +  &
                  REAL(MB(1)-MA1)*ALOGM2 + REAL(MA0)
          B2RDB = LOG(REAL(ABS(MB(2))+1)/REAL(ABS(MB2)+1))/0.69315 +  &
                  REAL(MB(1)-MB1)*ALOGM2 + REAL(MB0)
          MB2RD = NINT(MAX(0.0,MIN(B2RDA,B2RDB,  &
                  (NDIG-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)))
          IF (MB(2) == 0) THEN
              MB(0) = 0
          ELSE
              MB(0) = MIN(MAX(MA0,MB0),MB2RD)
          ENDIF
      ELSE
          MB(0) = MA0
      ENDIF

      JRSIGN = JRSSAV
      RETURN
      END SUBROUTINE FMADD2_R2

      SUBROUTINE FMADDI(MA,IVAL)

!  MA = MA + IVAL

!  Increment MA by one word integer IVAL.

!  This routine is faster than FMADD when IVAL is small enough so
!  that it can be added to a single word of MA without often causing
!  a carry.  Otherwise FMI2M and FMADD are used.

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK)
      INTEGER :: IVAL
      REAL (KIND(1.0D0)) :: MAEXP,MD2B,MKSUM
      INTEGER :: KPTMA

      NCALL = NCALL + 1
      IF (NTRACE /= 0) THEN
          NAMEST(NCALL) = 'FMADDI'
          CALL FMNTR(2,MA,MA,1,1)
          CALL FMNTRI(2,IVAL,0)
      ENDIF
      KFLAG = 0

      MAEXP = MA(1)
      IF (MAEXP <= 0 .OR. MAEXP > NDIG) GO TO 110
      KPTMA = INT(MAEXP) + 1
      IF (MA(-1) < 0) THEN
          MKSUM = MA(KPTMA) - IVAL
      ELSE
          MKSUM = MA(KPTMA) + IVAL
      ENDIF

      IF (MKSUM >= MBASE .OR. MKSUM < 0) GO TO 110
      IF (KPTMA == 2 .AND. MKSUM == 0) GO TO 110
      MA(KPTMA) = MKSUM
      GO TO 120

  110 CALL FMI2M(IVAL,M01)
      CALL FMADD2_R1(MA,M01)

  120 IF (KACCSW == 1) THEN
          MD2B = NINT((NDIG-1)*ALOGM2 + LOG(REAL(ABS(MA(2))+1))/0.69315)
          MA(0) = MIN(MA(0),MD2B)
      ENDIF
      IF (NTRACE /= 0) THEN
          CALL FMNTR(1,MA,MA,1,1)
      ENDIF
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMADDI

      SUBROUTINE FMADDN(MA,MB,NGUARD,NMWA)

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

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

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      INTEGER NGUARD,NMWA
      REAL (KIND(1.0D0)) :: MK,MR
      INTEGER J,K,KL,KP1,KP2,KPT,KSH,N1,N2,NK,NK1

      N1 = NDIG + 1

!             Check for an insignificant operand.

      MK = MA(1) - MB(1)
      IF (MK >= NDIG+2) THEN
          DO J = 1, N1
             MWA(J) = MA(J)
          ENDDO
          MWA(N1+1) = 0
          IF (KROUND == 0 .OR. (KROUND == 2 .AND. JRSIGN == -1) .OR.  &
             (KROUND == -1 .AND. JRSIGN == 1)) THEN
              MWA(N1) = MWA(N1) - 1
              GO TO 120
          ENDIF
          KFLAG = 1
          RETURN
      ENDIF
      K = INT(MK)
      IF (NGUARD <= 1) NMWA = N1 + 2

!             Subtract MB from MA.

      KP1 = MIN(N1,K+1)
      MWA(K+1) = 0
      DO J = 1, KP1
         MWA(J) = MA(J)
      ENDDO
      KP2 = K + 2

!             (Inner Loop)

      DO J = KP2, N1
         MWA(J) = MA(J) - MB(J-K)
      ENDDO

      N2 = NDIG + 2
      IF (N2-K <= 1) N2 = 2 + K
      NK = MIN(NMWA,N1+K)
      DO J = N2, NK
         MWA(J) = -MB(J-K)
      ENDDO
      NK1 = NK + 1
      DO J = NK1, NMWA
         MWA(J) = 0
      ENDDO

!             Normalize.  Fix the sign of any negative digit.

      IF (K > 0) THEN
          DO J = NMWA, KP2, -1
             IF (MWA(J) < 0) THEN
                 MWA(J) = MWA(J) + MBASE
                 MWA(J-1) = MWA(J-1) - 1
             ENDIF
          ENDDO

          KPT = KP2 - 1
  110     IF (MWA(KPT) < 0 .AND. KPT >= 3) THEN
              MWA(KPT) = MWA(KPT) + MBASE
              MWA(KPT-1) = MWA(KPT-1) - 1
              KPT = KPT - 1
              GO TO 110
          ENDIF
          GO TO 130
      ENDIF

  120 DO J = N1, 3, -1
         IF (MWA(J) < 0) THEN
             MWA(J) = MWA(J) + MBASE
             MWA(J-1) = MWA(J-1) - 1
         ENDIF
      ENDDO

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

  130 DO J = 2, NMWA
         IF (MWA(J) > 0) THEN
             KSH = J - 2
             GO TO 140
         ENDIF
      ENDDO
      MWA(1) = 0
      RETURN

  140 IF (KSH > 0) THEN
          KL = NMWA - KSH
          DO J = 2, KL
             MWA(J) = MWA(J+KSH)
          ENDDO
          DO J = KL+1, NMWA
             MWA(J) = 0
          ENDDO
          MWA(1) = MWA(1) - KSH
          IF (MK >= NDIG+2) THEN
              MWA(N1) = MBASE - 1
          ENDIF
      ENDIF

!             Round the result.

      MR = 2*MWA(NDIG+2) + 1
      IF (KROUND == -1 .OR. KROUND == 2) THEN
          CALL FMRND(MWA,NDIG,NGUARD,0)
      ELSE IF (MR >= MBASE) THEN
          IF (MR-1 > MBASE .AND. MWA(N1) < MBASE-1) THEN
              IF (KROUND /= 0 .OR. NCALL > 1) THEN
                  MWA(N1) = MWA(N1) + 1
                  MWA(N1+1) = 0
              ENDIF
          ELSE
              CALL FMRND(MWA,NDIG,NGUARD,0)
          ENDIF
      ENDIF

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

      IF (ABS(MA(1)-MB(1)) < NDIG) GO TO 150
      IF (ABS(MA(1)-MB(1)) > NDIG+1) THEN
          KFLAG = 1
          GO TO 150
      ENDIF

      N2 = NDIG + 4
      DO J = 3, N1
         IF (MWA(N2-J) /= MA(N2-J)) GO TO 150
      ENDDO
      IF (MWA(1) /= MA(1)) GO TO 150
      IF (MWA(2) /= ABS(MA(2))) GO TO 150
      KFLAG = 1

  150 RETURN
      END SUBROUTINE FMADDN

      SUBROUTINE FMADDP(MA,MB,NGUARD,NMWA)

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

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

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      INTEGER NGUARD,NMWA
      REAL (KIND(1.0D0)) :: MK,MKT,MR
      INTEGER J,K,KP,KP2,KPT,KSHIFT,N1,N2,NK

      N1 = NDIG + 1

!             Check for an insignificant operand.

      MK = MA(1) - MB(1)
      IF (MK >= NDIG+1) THEN
          MWA(1) = MA(1) + 1
          MWA(2) = 0
          DO J = 2, N1
             MWA(J+1) = MA(J)
          ENDDO
          MWA(N1+2) = 0
          IF ((KROUND ==  2 .AND. JRSIGN ==  1) .OR.  &
              (KROUND == -1 .AND. JRSIGN == -1)) THEN
              MWA(N1+2) = 1
              GO TO 120
          ENDIF
          KFLAG = 1
          RETURN
      ENDIF
      K = INT(MK)

!             Add MA and MB.

      MWA(1) = MA(1) + 1
      MWA(2) = 0
      DO J = 2, K+1
         MWA(J+1) = MA(J)
      ENDDO
      KP2 = K + 2

!             (Inner Loop)

      DO J = KP2, N1
         MWA(J+1) = MA(J) + MB(J-K)
      ENDDO
      N2 = NDIG + 2
      NK = MIN(NMWA,N1+K)
      DO J = N2, NK
         MWA(J+1) = MB(J-K)
      ENDDO
      DO J = NK+1, NMWA
         MWA(J+1) = 0
      ENDDO

!             Normalize.  Fix any digit not less than MBASE.

      IF (K == NDIG) GO TO 140

      IF (K > 0) THEN
          DO J = N1+1, KP2, -1
             IF (MWA(J) >= MBASE) THEN
                 MWA(J) = MWA(J) - MBASE
                 MWA(J-1) = MWA(J-1) + 1
             ENDIF
          ENDDO

          KPT = KP2 - 1
  110     IF (MWA(KPT) >= MBASE .AND. KPT >= 3) THEN
              MWA(KPT) = MWA(KPT) - MBASE
              MWA(KPT-1) = MWA(KPT-1) + 1
              KPT = KPT - 1
              GO TO 110
          ENDIF
          GO TO 120
      ENDIF

      DO J = N1+1, 3, -1
         IF (MWA(J) >= MBASE) THEN
             MWA(J) = MWA(J) - MBASE
             MWA(J-1) = MWA(J-1) + 1
         ENDIF
      ENDDO

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

  120 IF (MWA(2) >= MBASE) THEN
  130     KP = NMWA + 4
          DO J = 4, NMWA
             MWA(KP-J) = MWA(KP-J-1)
          ENDDO
          MKT = AINT (MWA(2)/MBASE)
          MWA(3) = MWA(2) - MKT*MBASE
          MWA(2) = MKT
          MWA(1) = MWA(1) + 1
          IF (MWA(2) >= MBASE) GO TO 130
      ENDIF

!             Round the result.

  140 KSHIFT = 0
      IF (MWA(2) == 0) KSHIFT = 1
      MR = 2*MWA(NDIG+2+KSHIFT) + 1
      IF (KROUND == -1 .OR. KROUND == 2) THEN
          CALL FMRND(MWA,NDIG,NGUARD,KSHIFT)
      ELSE IF (MR >= MBASE) THEN
          IF (MR-1 > MBASE .AND. MWA(N1+KSHIFT) < MBASE-1) THEN
              IF (KROUND /= 0 .OR. NCALL > 1) THEN
                  MWA(N1+KSHIFT) = MWA(N1+KSHIFT) + 1
                  MWA(N1+1+KSHIFT) = 0
              ENDIF
          ELSE
              CALL FMRND(MWA,NDIG,NGUARD,KSHIFT)
          ENDIF
      ENDIF

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

      IF (ABS(MA(1)-MB(1)) < NDIG) GO TO 150
      IF (KSHIFT == 0) GO TO 150
      IF (ABS(MA(1)-MB(1)) > NDIG+1) THEN
          KFLAG = 1
          GO TO 150
      ENDIF

      N2 = NDIG + 4
      DO J = 3, N1
         IF (MWA(N2-J+1) /= MA(N2-J)) GO TO 150
      ENDDO
      IF (MWA(1) /= MA(1)+1) GO TO 150
      IF (MWA(3) /= ABS(MA(2))) GO TO 150
      KFLAG = 1

  150 RETURN
      END SUBROUTINE FMADDP

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

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

!  KROUTN - Name of the subroutine that 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

      USE FMVALS
      IMPLICIT NONE
      CHARACTER(6) :: KROUTN
      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      INTEGER NARGS,KRESLT

      REAL (KIND(1.0D0)) :: MBS
      INTEGER J,KWRNSV,NCATMA,NCATMB,NDS

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

!             For example, row 7 column 2 of this array initialization
!             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.

      INTEGER :: KADD(15,15) = RESHAPE(  (/  &
                 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   /)  &
        , (/ 15,15 /) )

      INTEGER :: KMPY(15,15) = RESHAPE(  (/  &
                 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   /)  &
        , (/ 15,15 /) )

      INTEGER :: KDIV(15,15) = RESHAPE(  (/  &
                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   /)  &
        , (/ 15,15 /) )

      INTEGER :: KPWR(15,15) = RESHAPE(  (/  &
                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   /)  &
        , (/ 15,15 /) )

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

      KRESLT = 12
      KFLAG = -4
      IF (MA(1) == MUNKNO) RETURN
      IF (NARGS == 2) THEN
          IF (MB(1) == MUNKNO) RETURN
      ENDIF
      IF (MBLOGS /= MBASE) CALL FMCONS
      KFLAG = 0
      NAMEST(NCALL) = KROUTN

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

      IF (NCALL > 1 .AND. KDEBUG == 0) GO TO 130

!             Check NDIG.

      IF (NDIG < 2 .OR. NDIG > NDIGMX) THEN
          KFLAG = -1
          CALL FMWARN
          NDS = NDIG
          IF (NDIG < 2) NDIG = 2
          IF (NDIG > NDIGMX) NDIG = NDIGMX
          WRITE (KW,  &
                 "(' NDIG was',I10,'.  It has been changed to',I10,'.')"  &
                ) NDS,NDIG
          RETURN
      ENDIF

!             Check MBASE.

      IF (MBASE < 2 .OR. MBASE > MXBASE) THEN
          KFLAG = -2
          CALL FMWARN
          MBS = MBASE
          IF (MBASE < 2) MBASE = 2
          IF (MBASE > MXBASE) MBASE = MXBASE
          WRITE (KW,  &
                 "(' MBASE was',I10,'.  It has been changed to',I10,'.')"  &
                ) INT(MBS),INT(MBASE)
          CALL FMCONS
          RETURN
      ENDIF

!             Check exponent range.

      IF (MA(1) > MXEXP+1 .OR. MA(1) < -MXEXP) THEN
          IF (ABS(MA(1)) /= MEXPOV .OR. ABS(MA(2)) /= 1) THEN
              KFLAG = -3
              CALL FMWARN
              RETURN
          ENDIF
      ENDIF
      IF (NARGS == 2) THEN
          IF (MB(1) > MXEXP+1 .OR. MB(1) < -MXEXP) THEN
              IF (ABS(MB(1)) /= MEXPOV .OR. ABS(MB(2)) /= 1) THEN
                  KFLAG = -3
                  CALL FMWARN
                  RETURN
              ENDIF
          ENDIF
      ENDIF

!             Check for properly normalized digits in the
!             input arguments.

      IF (ABS(MA(1)-INT(MA(1))) /= 0) KFLAG = 1
      IF (MA(2) <= (-1) .OR. MA(2) >= MBASE .OR.  &
          ABS(MA(2)-INT(MA(2))) /= 0) KFLAG = 2
      IF (KDEBUG == 0) GO TO 110
      DO J = 3, NDIG+1
         IF (MA(J) < 0 .OR. MA(J) >= MBASE .OR.  &
             ABS(MA(J)-INT(MA(J))) /= 0) THEN
             KFLAG = J
             GO TO 110
         ENDIF
      ENDDO
  110 IF (KFLAG /= 0) THEN
          J = KFLAG
          KFLAG = -4
          KWRNSV = KWARN
          IF (KWARN >= 2) KWARN = 1
          CALL FMWARN
          KWARN = KWRNSV
          IF (KWARN >= 1) THEN
              WRITE (KW,*) ' First invalid array element:  MA(',  &
                           J,') = ',MA(J)
          ENDIF
          KFLAG = -4
          IF (KWARN >= 2) THEN
              STOP
          ENDIF
          RETURN
      ENDIF
      IF (NARGS == 2) THEN
          IF (ABS(MB(1)-INT(MB(1))) /= 0) KFLAG = 1
          IF (MB(2) <= (-1) .OR. MB(2) >= MBASE .OR.  &
              ABS(MB(2)-INT(MB(2))) /= 0) KFLAG = 2
          IF (KDEBUG == 0) GO TO 120
          DO J = 3, NDIG+1
             IF (MB(J) < 0 .OR. MB(J) >= MBASE .OR.  &
                 ABS(MB(J)-INT(MB(J))) /= 0) THEN
                 KFLAG = J
                 GO TO 120
             ENDIF
          ENDDO
  120     IF (KFLAG /= 0) THEN
              J = KFLAG
              KFLAG = -4
              KWRNSV = KWARN
              IF (KWARN >= 2) KWARN = 1
              CALL FMWARN
              KWARN = KWRNSV
              IF (KWARN >= 1) THEN
                  WRITE (KW,*) ' First invalid array element:  MB(',  &
                               J,') = ',MB(J)
              ENDIF
              KFLAG = -4
              IF (KWARN >= 2) THEN
                  STOP
              ENDIF
              RETURN
          ENDIF
      ENDIF

!             Check for special cases.

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

      IF (KROUTN == 'FMADD ') THEN
          KRESLT = KADD(NCATMB,NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMSUB ') THEN
          IF (NCATMB < 16) NCATMB = 16 - NCATMB
          KRESLT = KADD(NCATMB,NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMMPY ') THEN
          KRESLT = KMPY(NCATMB,NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMDIV ') THEN
          KRESLT = KDIV(NCATMB,NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMPWR ') THEN
          KRESLT = KPWR(NCATMB,NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMSQRT') THEN
          KRESLT = KSQRT(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMEXP ') THEN
          KRESLT = KEXP(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMLN  ') THEN
          KRESLT = KLN(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMSIN ') THEN
          KRESLT = KSIN(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMCOS ') THEN
          KRESLT = KCOS(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMTAN ') THEN
          KRESLT = KTAN(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMASIN') THEN
          KRESLT = KASIN(NCATMA)
          IF ((NCATMA == 7.OR.NCATMA == 9) .AND. KRAD == 0) KRESLT = 12
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMACOS') THEN
          KRESLT = KACOS(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMATAN') THEN
          KRESLT = KATAN(NCATMA)
          IF ((NCATMA == 7.OR.NCATMA == 9) .AND. KRAD == 0) KRESLT = 12
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMSINH') THEN
          KRESLT = KSINH(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMCOSH') THEN
          KRESLT = KCOSH(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMTANH') THEN
          KRESLT = KTANH(NCATMA)
          GO TO 140
      ENDIF

      IF (KROUTN == 'FMLG10') THEN
          KRESLT = KLG10(NCATMA)
          GO TO 140
      ENDIF

      KRESLT = 0
      RETURN

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

      SUBROUTINE FMASIN(MA,MB)

!  MB = ARCSIN(MA)

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      REAL (KIND(1.0D0)) :: MACCA,MACMAX,MXSAVE
      INTEGER J,K,KASAVE,KOVUN,KRESLT,NDSAVE

      IF (MBLOGS /= MBASE) CALL FMCONS
      IF (ABS(MA(1)) > MEXPAB .OR. MA(1) > 0 .OR. MA(2) == 0) THEN
          CALL FMENTR('FMASIN',MA,MA,1,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,  &
                      KOVUN)
          IF (KRESLT /= 0) RETURN
      ELSE
          NCALL = NCALL + 1
          NAMEST(NCALL) = 'FMASIN'
          IF (NTRACE /= 0) CALL FMNTR(2,MA,MA,1,1)
          KOVUN = 0
          IF (MA(1) == MEXPOV .OR. MA(1) == MEXPUN) KOVUN = 1
          NDSAVE = NDIG
          IF (NCALL == 1) THEN
              K = MAX(NGRD52-1,2)
              NDIG = MAX(NDIG+K,2)
              IF (NDIG > NDG2MX) THEN
                  KFLAG = -9
                  CALL FMWARN
                  NDIG = NDSAVE
                  KRESLT = 12
                  CALL FMRSLT(MA,MA,MB,KRESLT)
                  IF (NTRACE /= 0) CALL FMNTR(1,MB,MB,1,1)
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
          KASAVE = KACCSW
          KACCSW = 0
          MXSAVE = MXEXP
          MXEXP = MXEXP2
      ENDIF

      MACCA = MA(0)
      CALL FMEQ2(MA,MB,NDSAVE,NDIG)
      MB(0) = NINT(NDIG*ALOGM2)

!             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_R2(M03,M04)
      CALL FMSQRT_R1(M04)
      CALL FMDIV_R1(MB,M04)

      CALL FMATAN(MB,M13)
      CALL FMEQ(M13,MB)

!             Round the result and return.

      MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)
      MB(0) = MIN(MB(0),MACCA,MACMAX)
      DO J = -1, NDIG+1
         M01(J) = MB(J)
      ENDDO
      CALL FMEXIT(M01,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
      RETURN
      END SUBROUTINE FMASIN

      SUBROUTINE FMATAN(MA,MB)

!  MB = ARCTAN(MA)

      USE FMVALS
      IMPLICIT NONE

      DOUBLE PRECISION X,XM
      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK)
      INTEGER NSTACK(19)
      REAL (KIND(1.0D0)) :: MA1,MACCA,MACMAX,MAS,MXSAVE
      INTEGER J,K,KASAVE,KOVUN,KRESLT,KRSAVE,KST,KWRNSV,NDSAV1,NDSAVE,  &
              NDSV

      IF (MBLOGS /= MBASE) CALL FMCONS
      IF (ABS(MA(1)) > MEXPAB .OR. MA(2) == 0) THEN
          CALL FMENTR('FMATAN',MA,MA,1,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,  &
                      KOVUN)
          IF (KRESLT /= 0) RETURN
      ELSE
          NCALL = NCALL + 1
          NAMEST(NCALL) = 'FMATAN'
          IF (NTRACE /= 0) CALL FMNTR(2,MA,MA,1,1)
          KOVUN = 0
          IF (MA(1) == MEXPOV .OR. MA(1) == MEXPUN) KOVUN = 1
          NDSAVE = NDIG
          IF (NCALL == 1) THEN
              K = MAX(NGRD52-1,2)
              NDIG = MAX(NDIG+K,2)
              IF (NDIG > NDG2MX) THEN
                  KFLAG = -9
                  CALL FMWARN
                  NDIG = NDSAVE
                  KRESLT = 12
                  CALL FMRSLT(MA,MA,MB,KRESLT)
                  IF (NTRACE /= 0) CALL FMNTR(1,MB,MB,1,1)
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
          KASAVE = KACCSW
          KACCSW = 0
          MXSAVE = MXEXP
          MXEXP = MXEXP2
      ENDIF

      MACCA = MA(0)
      CALL FMEQ2(MA,M05,NDSAVE,NDIG)
      M05(0) = NINT(NDIG*ALOGM2)

!             If MA >= 1 work with 1/MA.

      MA1 = MA(1)
      MAS = MA(-1)
      M05(-1) = 1
      IF (MA1 >= 1) THEN
          CALL FMI2M(1,MB)
          CALL FMDIV_R2(MB,M05)
      ENDIF

      KRSAVE = KRAD
      KRAD = 1
      KWRNSV = KWARN

      X = M05(1)
      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 >= 1 .OR. KRSAVE == 0) THEN
          IF (MBSPI /= MBASE .OR. NDIGPI < NDIG)  THEN
              NDSV = NDIG
              NDIG = MIN(NDIG+2,NDG2MX)
              NCALL = NCALL + 1
              NAMEST(NCALL) = 'NOEQ  '
              CALL FMPI(MPISAV)
              NCALL = NCALL - 1
              NDIG = NDSV
          ENDIF
      ENDIF

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

      IF (X*DLOGMB < -5.0D0*LOG(XM)) THEN
          KWARN = 0
          CALL FMEQ(M05,MB)
          IF (MB(1) <= -NDIG) GO TO 120
          CALL FMSQR(M05,M06)
          J = 3
          NDSAV1 = NDIG

  110     CALL FMMPY_R1(M05,M06)
          IF (M05(1) /= MUNKNO .AND. M05(2) /= 0) M05(-1) = -M05(-1)
          CALL FMDIVI(M05,J,M03)
          NDIG = NDSAV1
          CALL FMADD_R1(MB,M03)
          IF (KFLAG /= 0) THEN
              KFLAG = 0
              GO TO 120
          ENDIF
          NDIG = NDSAV1 - INT((MB(1)-M03(1)))
          IF (NDIG < 2) NDIG = 2
          J = J + 2
          GO TO 110
      ELSE

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

!             Newton iteration.

          DO J = 1, KST
             NDIG = NSTACK(J)
             CALL FMSIN(MB,M06)
             CALL FMSQR(M06,M03)
             CALL FMI2M(1,M04)
             CALL FMSUB_R2(M04,M03)
             CALL FMSQRT(M03,M04)
             CALL FMDIV_R2(M06,M04)
             CALL FMSUB_R1(M04,M05)
             CALL FMMPY_R2(M03,M04)
             CALL FMSUB_R1(MB,M04)
          ENDDO
          MB(0) = NINT(NDIG*ALOGM2)
      ENDIF

!             If MA >= 1 use pi/2 - ATAN(1/MA)

  120 IF (MA1 >= 1) THEN
          CALL FMDIVI(MPISAV,2,M06)
          CALL FMSUB_R2(M06,MB)
      ENDIF

!             Convert to degrees if necessary, round and return.

      KRAD = KRSAVE
      IF (KRAD == 0) THEN
          CALL FMMPYI_R1(MB,180)
          CALL FMDIV_R1(MB,MPISAV)
      ENDIF
      IF (MB(1) /= MUNKNO .AND. MB(2) /= 0 .AND. MAS < 0) MB(-1) = -MB(-1)

      IF (KFLAG == 1) KFLAG = 0
      KWARN = KWRNSV
      MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)
      MB(0) = MIN(MB(0),MACCA,MACMAX)
      DO J = -1, NDIG+1
         M01(J) = MB(J)
      ENDDO
      CALL FMEXIT(M01,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
      RETURN
      END SUBROUTINE FMATAN

      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.

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK),MC(-1:LUNPCK)
      REAL (KIND(1.0D0)) :: MACCA,MACCB,MACMAX,MXEXP1,MXSAVE
      INTEGER J,JQUAD,K,KASAVE,KOVUN,KRESLT,KWRNSV,NDSAVE

      IF (MBLOGS /= MBASE) CALL FMCONS
      IF (ABS(MA(1)) > MEXPAB .OR. ABS(MB(1)) > MEXPAB) THEN
          CALL FMENTR('FMATN2',MA,MB,2,1,MC,KRESLT,NDSAVE,MXSAVE,KASAVE,  &
                      KOVUN)
          IF (KRESLT /= 0) RETURN
      ELSE
          NCALL = NCALL + 1
          NAMEST(NCALL) = 'FMATN2'
          IF (NTRACE /= 0) CALL FMNTR(2,MA,MB,2,1)
          KOVUN = 0
          IF (MA(1) == MEXPOV .OR. MA(1) == MEXPUN) KOVUN = 1
          NDSAVE = NDIG
          IF (NCALL == 1) THEN
              K = MAX(NGRD52-1,2)
              NDIG = MAX(NDIG+K,2)
              IF (NDIG > NDG2MX) THEN
                  KFLAG = -9
                  CALL FMWARN
                  NDIG = NDSAVE
                  KRESLT = 12
                  CALL FMRSLT(MA,MB,MC,KRESLT)
                  IF (NTRACE /= 0) CALL FMNTR(1,MC,MC,1,1)
                  NCALL = NCALL - 1
                  RETURN
              ENDIF
          ENDIF
          KASAVE = KACCSW
          KACCSW = 0
          MXSAVE = MXEXP
          MXEXP = MXEXP2
      ENDIF

      KWRNSV = KWARN
      KWARN = 0

      MACCA = MA(0)
      MACCB = MB(0)
      CALL FMEQ2(MA,M01,NDSAVE,NDIG)
      M01(0) = NINT(NDIG*ALOGM2)
      CALL FMEQ2(MB,M02,NDSAVE,NDIG)
      M02(0) = NINT(NDIG*ALOGM2)

!             Check for special cases.

      IF (MA(1) == MUNKNO .OR. MB(1) == MUNKNO .OR.  &
         (MA(2) == 0 .AND. MB(2) == 0)) THEN
          CALL FMST2M('UNKNOWN',MC)
          KFLAG = -4
          GO TO 110
      ENDIF

      IF (MB(2) == 0 .AND. MA(-1) > 0) THEN
          IF (KRAD == 0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI_R1(MC,2)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(2) == 0 .AND. MA(-1) < 0) THEN
          IF (KRAD == 0) THEN
              CALL FMI2M(-90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI_R1(MC,-2)
          ENDIF
          GO TO 110
      ENDIF

      MXEXP1 = INT(MXEXP2/2.01D0)
      IF (MA(1) == MEXPOV .AND. MB(1) < MXEXP1-NDIG-2) THEN
          IF (KRAD == 0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI_R1(MC,2)
          ENDIF
          IF (M01(-1) < 0) MC(-1) = -1
          GO TO 110
      ENDIF

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

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

      IF (MB(1) == MEXPUN .AND. MA(2) == 0) THEN
          IF (MB(-1) < 0) THEN
              IF (KRAD == 0) THEN
                  CALL FMI2M(180,MC)
              ELSE
                  CALL FMPI(MC)
              ENDIF
          ELSE
              CALL FMI2M(0,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(1) == MEXPUN .AND. (-MA(1)) < MXEXP1-NDIG-2) THEN
          IF (KRAD == 0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI_R1(MC,2)
          ENDIF
          IF (M01(-1) < 0) MC(-1) = -1
          GO TO 110
      ENDIF

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

      IF (MA(-1) >= 0 .AND. MB(-1) > 0) JQUAD = 1
      IF (MA(-1) >= 0 .AND. MB(-1) < 0) JQUAD = 2
      IF (MA(-1) < 0 .AND. MB(-1) < 0) JQUAD = 3
      IF (MA(-1) < 0 .AND. MB(-1) > 0) JQUAD = 4

      CALL FMDIV(M01,M02,MC)
      MC(-1) = 1
      CALL FMATAN(MC,M13)
      CALL FMEQ(M13,MC)

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

      IF ((JQUAD == 3 .OR. JQUAD == 4) .AND. MC(1) /= MUNKNO .AND.  &
          MC(2) /= 0) MC(-1) = -MC(-1)

!             Round the result and return.

  110 IF (KFLAG == 1) KFLAG = 0
      KWARN = KWRNSV
      MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315)
      MC(0) = MIN(MC(0),MACCA,MACCB,MACMAX)
      DO J = -1, NDIG+1
         M01(J) = MC(J)
      ENDDO
      CALL FMEXIT(M01,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
      RETURN
      END SUBROUTINE FMATN2

      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.

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK)

      INTEGER J,N1

      NCALL = NCALL + 1
      NAMEST(NCALL) = 'FMBIG '

      IF (MBLOGS /= MBASE) CALL FMCONS
      KFLAG = 0
      N1 = NDIG + 1
      DO J = 2, N1
         MA(J) = MBASE - 1
      ENDDO
      MA(1) = MXEXP + 1
      MA(0) = NINT(NDIG*ALOGM2)
      MA(-1) = 1

      IF (NTRACE /= 0) CALL FMNTR(1,MA,MA,1,1)
      NCALL = NCALL - 1
      RETURN
      END SUBROUTINE FMBIG

      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) >= 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) >= MAXEXP-NDIG+1 )
!   8.          0
!   9.         +UN                UN stands for underflowed results.
!  10.   (+UN   , +UNTH)             ( MA(1) <= -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) <= -MAXEXP+NDIG-1 )
!  16.       UNKNOWN

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK)
      INTEGER NCAT

      REAL (KIND(1.0D0)) :: MA2,MXEXP1
      INTEGER J,NLAST

!             Check for special symbols.

      NCAT = 16
      IF (MA(1) == MUNKNO) RETURN

      IF (MA(1) == MEXPOV) THEN
          NCAT = 15
          IF (MA(-1) < 0) NCAT = 1
          RETURN
      ENDIF

      IF (MA(1) == MEXPUN) THEN
          NCAT = 9
          IF (MA(-1) < 0) NCAT = 7
          RETURN
      ENDIF

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

!             Check for +1 or -1.

      MA2 = ABS(MA(2))
      IF (MA(1) == 1 .AND. MA2 == 1) THEN
          NLAST = NDIG + 1
          IF (NLAST >= 3) THEN
              DO J = 3, NLAST
                 IF (MA(J) /= 0) GO TO 110
              ENDDO
          ENDIF
          NCAT = 12
          IF (MA(-1) < 0) NCAT = 4
          RETURN
      ENDIF

  110 MXEXP1 = INT(MXEXP)
      IF (MA(1) >= MXEXP1-NDIG+2) THEN
          NCAT = 14
          IF (MA(-1) < 0) NCAT = 2
          RETURN
      ENDIF

      IF (MA(1) >= 1) THEN
          NCAT = 13
          IF (MA(-1) < 0) NCAT = 3
          RETURN
      ENDIF

      IF (MA(1) >= -MXEXP1+NDIG) THEN
          NCAT = 11
          IF (MA(-1) < 0) NCAT = 5
          RETURN
      ENDIF

      IF (MA(1) >= -MXEXP1) THEN
          NCAT = 10
          IF (MA(-1) < 0) NCAT = 6
          RETURN
      ENDIF

      RETURN
      END SUBROUTINE FMCAT

      SUBROUTINE FMCHSH(MA,MB,MC)

!  MB = COSH(MA),    MC = SINH(MA)

!  If both the hyperbolic sine and cosine are needed, this routine
!  is faster than calling both FMCOSH and FMSINH.

!  MB and MC must be distinct arrays.

      USE FMVALS
      IMPLICIT NONE

      REAL (KIND(1.0D0)) :: MA(-1:LUNPCK),MB(-1:LUNPCK),MC(-1:LUNPCK)
      REAL (KIND(1.0D0)) :: MACCA,MACMAX,MAS,MXSAVE
      INTEGER J,K,KASAVE,KOVUN,KRESLT,KWRNSV,NCSAVE,NDSAVE

      IF (MBLOGS /= MBASE) CALL FMCONS
      MACCA = MA(0)
      MAS = MA(-1)
      IF (ABS(MA(1)) > MEXPAB) THEN
          NCSAVE = NCALL
          CALL FMENTR('FMCHSH',MA,MA,1,1,MB,KRESLT,NDSAVE,MXSAVE,KASAVE,  &
                      KOVUN)
          IF (MA(1) == MUNKNO) KOVUN = 2
          NCALL = NCSAVE + 1
          CALL FMEQ2(MA,M04,NDSAVE,NDIG)
          M04(0) = NINT(NDIG*ALOGM2)
          M04(-1) = 1
          CALL FMCOSH(M04,MB)
          CALL FMSINH(M04,MC)
          GO TO 110
      ELSE
          NCALL = NCALL + 1
          NAMEST(NCALL) = 'FMCHSH'
          IF (NTRACE /= 0) CALL FMNTR(2,MA,MA,1,1)
          KOVUN = 0
          IF (MA(1) == MEXPOV .OR. MA(1) == MEXPUN) KOVUN = 1
          NDSAVE = NDIG
          IF (NCALL == 1) THEN
              K = MAX(NGRD52,2)
              NDIG = MAX(NDIG+K,2)
              IF (NDIG > NDG2MX) THEN
                  NCALL = NCALL - 1
                  NDIG = NDSAVE
                  CALL FMEQ(MA,M04)
                  CALL FMCOSH(M04,MB)
                  CALL FMSINH(M04,MC)
                  KFLAG = -9
                  RETURN
              ENDIF
          ENDIF
          KASAVE = KACCSW
          KACCSW = 0
          MXSAVE = MXEXP
          MXEXP = MXEXP2
      ENDIF

      CALL FMEQ2(MA,M04,NDSAVE,NDIG)
      M04(0) = NINT(NDIG*ALOGM2)
      M04(-1) = 1

      K = 1
      IF (M04(1) == 0 .AND. M04(2) /= 0) THEN
          IF (MBASE/M04(2) >= 100) K = 2
      ENDIF
      IF (M04(1) >= 0 .AND. M04(2) /= 0 .AND. K == 1) THEN
          CALL FMCOSH(M04,MB)
          IF (MB(1) > NDIG) THEN
              CALL FMEQ(MB,MC)
              GO TO 110
          ENDIF
          CALL FMSQR(MB,M03)
          CALL FMI2M(-1,M02)
          CALL FMADD_R1(M03,M02)
          CALL FMSQRT(M03,MC)
      ELSE
          CALL FMSINH(M04,MC)
          CALL FMSQR(MC,M03)
          CALL FMI2M(1,M02)
          CALL FMADD_R1(M03,M02)
          CALL FMSQRT(M03,MB)
      ENDIF

!             Round and return.

  110 MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MC(2))+1))/0.69315)
      MC(0) = MIN(MC(0),MACCA,MACMAX)
      IF (MAS < 0 .AND. MC(1) /= MUNKNO .AND. MC(2) /= 0) MC(-1) = -MC(-1)
      CALL FMEQ2_R1(MC,NDIG,NDSAVE)
      MACMAX = NINT((NDSAVE-1)*ALOGM2 + LOG(REAL(ABS(MB(2))+1))/0.69315)
      MB(0) = MIN(MB(0),MACCA,MACMAX)
      IF (KOVUN == 2) THEN
          KWRNSV = KWARN
          KWARN = 0
      ENDIF
      DO J = -1, NDIG+1
         M01(J) = MB(J)
      ENDDO
      CALL FMEXIT(M01,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
      IF (KOVUN == 2) THEN
          KWARN = KWRNSV
      ENDIF
      IF (NTRACE /= 0) THEN
          IF (ABS(NTRACE) >= 1 .AND. NCALL+1 <= LVLTRC) THEN
              IF (NTRACE < 0) THEN
                  CALL FMNTRJ(MC,NDIG)
              ELSE
                  CALL FMPRNT(MC)