xp_module.f90

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

   CALL FMMPY(MA%DIGITS,MB%DIGITS,XP_MPY%DIGITS) ! XP_MPY = MA*MB
END FUNCTION

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

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

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

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

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

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

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

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

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

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

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

SUBROUTINE XP_SET (NPREC)

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

END SUBROUTINE

SUBROUTINE FMABS(MA,MB)

!  MB = ABS(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

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

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

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

SUBROUTINE FMACOS(MA,MB)

!  MB = ACOS(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

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

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

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

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

      CALL FMATAN(MB,MB)

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

!             Round the result and return.

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

SUBROUTINE FMADD(MA,MB,MC)

!  MC = MA + MB

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

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

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

      CALL FMADD2(MA,MB,MC)

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

SUBROUTINE FMADD2(MA,MB,MC)

!  Internal addition routine.  MC = MA + MB

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

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

      KFLAG = 0
      N1 = NDIG + 1

!             NGUARD is the number of guard digits used.

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

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

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

!             See which one is larger in absolute value.

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

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

      JCOMP = 2

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


!             Round the result.

      CALL FMRND(NDIG,NGUARD,0)

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

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

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

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

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

END SUBROUTINE

SUBROUTINE FMADDN(MA,MB,NGUARD,NMWA)

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

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

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      N1 = NDIG + 1

!             Check for an insignificant operand.

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

!             Subtract MB from MA.

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

!             Normalize.  Fix the sign of any negative digit.

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

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

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

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

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

END SUBROUTINE

SUBROUTINE FMADDP(MA,MB,NMWA)

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

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

      DIMENSION MA(LUNPCK),MB(LUNPCK)

      N1 = NDIG + 1

!             Check for an insignificant operand.

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

!             Add MA and MB.

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