xp_module.f90

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

      ENDIF
      NK1 = NK + 1
      IF (NK1.LE.NMWA) THEN
          DO 150 J = NK1, NMWA
             MWA(J) = 0
  150     CONTINUE
      ENDIF

!             Normalize.  Fix any digit not less than JBASE.

      IF (K.EQ.NDIG) RETURN

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

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

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

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

END SUBROUTINE

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

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

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

!  Result codes:

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

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

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

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

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

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

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

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

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

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

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

      IF (NCALL.GT.1) GO TO 130

!             Check NDIG.

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

!             Check JBASE.

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

!             Check exponent range.

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

!             Check for special cases.

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

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

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

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

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

      CASE (39)
          KRESLT = KSQRT(NCATMA)

      CASE (15)
          KRESLT = KEXP(NCATMA)

      CASE (21)
          KRESLT = KLN(NCATMA)

      CASE (36)
          KRESLT = KSIN(NCATMA)

      CASE (9)
          KRESLT = KCOS(NCATMA)

      CASE (41)
          KRESLT = KTAN(NCATMA)

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

      CASE (2)
          KRESLT = KACOS(NCATMA)

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

      CASE (37)
          KRESLT = KSINH(NCATMA)

      CASE (10)
          KRESLT = KCOSH(NCATMA)

      CASE (42)
          KRESLT = KTANH(NCATMA)

      CASE (20)
          GO TO 140

      CASE DEFAULT
      KRESLT = 0

      END SELECT

      RETURN

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

SUBROUTINE FMASIN(MA,MB)

!  MB = ARCSIN(MA)

      DIMENSION MA(LUNPCK),MB(LUNPCK)

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

      CALL FMEQU(MA,MB,NDSAVE,NDIG)

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

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

      CALL FMATAN(MB,MB)

!             Round the result and return.

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

SUBROUTINE FMATAN(MA,MB)

!  MB = ARCTAN(MA)

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

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

      CALL FMEQU(MA,M05,NDSAVE,NDIG)

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

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

      KRSAVE = KRAD
      KRAD = 1
      KWSV = KWARN

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

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

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

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

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

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

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

!             Newton iteration.

          DO 120 J = 1, KST