xp_module.f90

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

!                    compute correctly rounded results.
!                    In the definition of LUNPCK The '6/5' condition
!                    allows for converting from a large base to the
!                    (smaller) largest power of ten base for output
!                    conversion.
!                    The '+ 20' condition allows for the need to carry
!                    many guard digits when using a small base like 2.

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

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

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

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

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

!             KFLAG is the flag for error conditions.

   INTEGER :: KFLAG = 0

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

   INTEGER :: NTRACE = 0

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

  INTEGER ::  LVLTRC = 1

!             NCALL is the call stack pointer.

  INTEGER, PRIVATE ::  NCALL = 0

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

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

  INTEGER, PRIVATE :: NDIGPI = 0

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

  INTEGER, PRIVATE :: NJBPI = 0

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

  INTEGER, PRIVATE :: NDIGE = 0

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

  INTEGER, PRIVATE :: NJBE = 0

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

  INTEGER, PRIVATE :: NDIGLB = 0

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

  INTEGER, PRIVATE :: NJBLB = 0

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

  INTEGER, PRIVATE :: NDIGLI = 0

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

  INTEGER, PRIVATE :: NJBLI = 0

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

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

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

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

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

  INTEGER, PRIVATE :: KARGSW = 0

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

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

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

  INTEGER, PRIVATE, PARAMETER :: KEXPOV = -KEXPUN

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

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

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

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

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

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

!             JFORM1 indicates the format used by FMOUT.

   INTEGER :: JFORM1 = 1

!             JFORM2 indicates the number of digits used in FMOUT.

   INTEGER :: JFORM2 = INITIAL_NPREC

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

   INTEGER :: KRAD = 1

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

   INTEGER :: KWARN = 1

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

   INTEGER :: KROUND = 1

CONTAINS

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

FUNCTION XP_MPY (MA, MB)
   TYPE (XP_REAL), INTENT (IN) :: MA, MB
   TYPE (XP_REAL) :: XP_MPY