xp_module.f90

任意多精度计算模块xp_module 很好用的

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

MODULE XP_MODULE

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

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

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

!     FM 1.0                  David M Smith                  1-06-90

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

!  The XP_REAL package performs floating point extended precision arithmetic.


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

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

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

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

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

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

!                        1    -31415  92653590

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


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

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

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


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


!  NTRACE and LVLTRC control trace printout from the package.

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


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

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

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

!    KFLAG =  0     Normal operation.

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

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

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

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

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


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


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


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

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

   INTERFACE ABS
      MODULE PROCEDURE XP_ABS
   END INTERFACE

   INTERFACE ACOS
      MODULE PROCEDURE XP_ACOS
   END INTERFACE

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

   INTERFACE AINT
      MODULE PROCEDURE XP_AINT
   END INTERFACE

   INTERFACE ASIN
      MODULE PROCEDURE XP_ASIN
   END INTERFACE

   INTERFACE ATAN
      MODULE PROCEDURE XP_ATAN
   END INTERFACE

   INTERFACE ATAN2
      MODULE PROCEDURE XP_ATAN2
   END INTERFACE

   INTERFACE HUGE
      MODULE PROCEDURE XP_HUGE
   END INTERFACE

   INTERFACE OPERATOR (==)
      MODULE PROCEDURE XP_EQ
   END INTERFACE

   INTERFACE OPERATOR (/=)
      MODULE PROCEDURE XP_NE
   END INTERFACE

   INTERFACE OPERATOR (<=)
      MODULE PROCEDURE XP_LE
   END INTERFACE

   INTERFACE OPERATOR (>=)
      MODULE PROCEDURE XP_GE
   END INTERFACE

   INTERFACE OPERATOR (<)
      MODULE PROCEDURE XP_LT
   END INTERFACE

   INTERFACE OPERATOR (>)
      MODULE PROCEDURE XP_GT
   END INTERFACE

   INTERFACE COS
      MODULE PROCEDURE XP_COS
   END INTERFACE

   INTERFACE COSH
      MODULE PROCEDURE XP_COSH
   END INTERFACE

   INTERFACE DIM
      MODULE PROCEDURE XP_DIM
   END INTERFACE

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

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

   INTERFACE EXP
      MODULE PROCEDURE XP_EXP
   END INTERFACE

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

   INTERFACE LOG10
      MODULE PROCEDURE XP_LOG10
   END INTERFACE

   INTERFACE LOG
      MODULE PROCEDURE XP_LOG
   END INTERFACE

   INTERFACE DBLE
      MODULE PROCEDURE XP_M2DP
   END INTERFACE

   INTERFACE INT
      MODULE PROCEDURE XP_M2I
   END INTERFACE

   INTERFACE REAL
      MODULE PROCEDURE XP_M2SP
   END INTERFACE

   INTERFACE MAX
      MODULE PROCEDURE XP_MAX
   END INTERFACE

   INTERFACE MIN
      MODULE PROCEDURE XP_MIN
   END INTERFACE

   INTERFACE MOD
      MODULE PROCEDURE XP_MOD
   END INTERFACE

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

   INTERFACE ANINT
      MODULE PROCEDURE XP_ANINT
   END INTERFACE

   INTERFACE PI
      MODULE PROCEDURE XP_PI
   END INTERFACE

   INTERFACE SIGN
      MODULE PROCEDURE XP_SIGN
   END INTERFACE

   INTERFACE SIN
      MODULE PROCEDURE XP_SIN
   END INTERFACE

   INTERFACE SINH
      MODULE PROCEDURE XP_SINH
   END INTERFACE

   INTERFACE SQRT
      MODULE PROCEDURE XP_SQRT
   END INTERFACE

   INTERFACE OPERATOR (-)
      MODULE PROCEDURE XP_SUB
   END INTERFACE

   INTERFACE TAN
      MODULE PROCEDURE XP_TAN
   END INTERFACE

   INTERFACE TANH
      MODULE PROCEDURE XP_TANH
   END INTERFACE

!   INTERFACE PRINT
!      MODULE PROCEDURE XP_PRINT
!   END INTERFACE

!   PRIVATE XP_PRINT


!  FM numbers in unpacked format have size LUNPCK.


!  PORTABILITY NOTES.

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

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

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

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

      PARAMETER ( MXNDIG=256 , NBITS=32 ,  &

!             Define the array sizes:

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

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

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


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

!  The following contain information about saved constants.

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

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

      INTEGER, PRIVATE :: MJSUMS(LJSUMS)

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

      CHARACTER, PRIVATE :: MBUFF(LMBUFF)

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

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

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

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

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

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

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

   INTEGER :: KW = 6

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

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

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

!             MXNDG2 is the maximum value for NDIG which can be used
!                    internally.  Several of the FM routines may raise
!                    NDIG above MXNDIG temporarily, in order to