xp_module.f90

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

             NDIG = NSTACK(J)
             CALL FMSIN(MB,M06)
             CALL FMMPY(M06,M06,M03)
             CALL FMI2M(1,M04)
             CALL FMSUB(M04,M03,M03)
             CALL FMSQRT(M03,M04)
             CALL FMDIV(M06,M04,M04)
             CALL FMSUB(M04,M05,M04)
             CALL FMMPY(M03,M04,M04)
             CALL FMSUB(MB,M04,MB)
  120     CONTINUE
      ENDIF

!             If MA.GE.1 use pi/2 - ATAN(1/MA)

  130 IF (MA1.GE.1) THEN
          CALL FMDIVI(MPISAV,2,M06)
          CALL FMSUB(M06,MB,MB)
      ENDIF

!             Convert to degrees if necessary, round and return.

      KRAD = KRSAVE
      IF (KRAD.EQ.0) THEN
          CALL FMMPYI(MB,180,MB)
          CALL FMDIV(MB,MPISAV,MB)
      ENDIF
      IF (MA2.LT.0) MB(2) = -MB(2)

      IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(MB,MB,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

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.

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

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

      KARGSW = 0
      KWSV = KWARN
      KWARN = 0
      CALL FMEQU(MA,M01,NDSAVE,NDIG)
      CALL FMEQU(MB,M02,NDSAVE,NDIG)

!             Check for special cases.

      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO .OR.  &
         (MA(2).EQ.0 .AND. MB(2).EQ.0)) THEN
          CALL FMIM(0,MC)
          MC(1) = KUNKNO
          MC(2) = 1
          KFLAG = -4
          GO TO 110
      ENDIF

      IF (MB(2).EQ.0 .AND. MA(2).GT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(2).EQ.0 .AND. MA(2).LT.0) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(-90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,-2,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MA(1).EQ.KEXPOV .AND. MB(1).LT.MXSAVE-NDIG-2) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

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

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

      IF (MB(1).EQ.KEXPUN .AND. MA(2).EQ.0) THEN
          IF (MB(2).LT.0) THEN
              IF (KRAD.EQ.0) THEN
                  CALL FMI2M(180,MC)
              ELSE
                  CALL FMPI(MC)
              ENDIF
          ELSE
              CALL FMI2M(0,MC)
          ENDIF
          GO TO 110
      ENDIF

      IF (MB(1).EQ.KEXPUN .AND. (-MA(1)).LT.MXSAVE-NDIG-2) THEN
          IF (KRAD.EQ.0) THEN
              CALL FMI2M(90,MC)
          ELSE
              CALL FMPI(MC)
              CALL FMDIVI(MC,2,MC)
          ENDIF
          IF (M01(2).LT.0) MC(2) = -MC(2)
          GO TO 110
      ENDIF

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

      IF (MA(2).GE.0 .AND. MB(2).GT.0) JQUAD = 1
      IF (MA(2).GE.0 .AND. MB(2).LT.0) JQUAD = 2
      IF (MA(2).LT.0 .AND. MB(2).LT.0) JQUAD = 3
      IF (MA(2).LT.0 .AND. MB(2).GT.0) JQUAD = 4

      CALL FMDIV(M01,M02,MC)
      MC(2) = ABS(MC(2))
      CALL FMATAN(MC,MC)

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

      IF ((JQUAD.EQ.3 .OR. JQUAD.EQ.4) .AND. MC(1).NE.KUNKNO)  &
          MC(2) = -MC(2)

!             Round the result and return.

  110 IF (KFLAG.EQ.1) KFLAG = 0
      KWARN = KWSV
      CALL FMEXIT(MC,MC,NDSAVE,MXSAVE,KASAVE,KOVUN)
END SUBROUTINE

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.

      DIMENSION MA(LUNPCK)

      NCALL = NCALL + 1
      KSTACK(NCALL) = 7

      KFLAG = 0
      N1 = NDIG + 1
      DO 110 J = 2, N1
         MA(J) = JBASE - 1
  110 CONTINUE
      MA(1) = MXEXP + 1

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

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

      DIMENSION MA(LUNPCK)

!             Check for special symbols.

      NCAT = 16
      IF (MA(1).EQ.KUNKNO) RETURN

      IF (MA(1).EQ.KEXPOV) THEN
          NCAT = 15
          IF (MA(2).LT.0) NCAT = 1
          RETURN
      ENDIF

      IF (MA(1).EQ.KEXPUN) THEN
          NCAT = 9
          IF (MA(2).LT.0) NCAT = 7
          RETURN
      ENDIF

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

!             Check for +1 or -1.

      MA2 = ABS(MA(2))
      IF (MA(1).EQ.1 .AND. MA2.EQ.1) THEN
          NLAST = NDIG + 1
          IF (NLAST.GE.3) THEN
              DO 110 J = 3, NLAST
                 IF (MA(J).NE.0) GO TO 120
  110         CONTINUE
          ENDIF
          NCAT = 12
          IF (MA(2).LT.0) NCAT = 4
          RETURN
      ENDIF

  120 IF (MA(1).GE.MXEXP-NDIG+1) THEN
          NCAT = 14
          IF (MA(2).LT.0) NCAT = 2
          RETURN
      ENDIF

      IF (MA(1).GE.1) THEN
          NCAT = 13
          IF (MA(2).LT.0) NCAT = 3
          RETURN
      ENDIF

      IF (MA(1).GE.-MXEXP+NDIG) THEN
          NCAT = 11
          IF (MA(2).LT.0) NCAT = 5
          RETURN
      ENDIF

      IF (MA(1).GE.-MXEXP) THEN
          NCAT = 10
          IF (MA(2).LT.0) NCAT = 6
          RETURN
      ENDIF

END SUBROUTINE

FUNCTION FMCOMP(MA,LREL,MB)
   LOGICAL FMCOMP
   DIMENSION MA(LUNPCK),MB(LUNPCK)
   CHARACTER *2 LREL

!  Logical comparison of FM numbers MA and MB.

!  LREL is a CHARACTER *2 description of the comparison to be done:
!  LREL = 'EQ' returns FMCOMP = .TRUE. if MA.EQ.MB
!       = 'NE', 'GE', 'GT', 'LE', 'LT' also work like a logical IF.

!  For comparisons involving 'UNKNOWN' or two identical special symbols
!  such as +OVERFLOW,'EQ',+OVERFLOW FMCOMP is returned FALSE and a
!  KFLAG = -4 error condition is returned.

      NCALL = NCALL + 1
      KSTACK(NCALL) = 8

      IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.2) THEN
          WRITE (KW,110)
  110     FORMAT(' INPUT TO FMCOMP')

          IF (NTRACE.GT.0) THEN
              CALL FMPRNT(MA)
              WRITE (KW,120) LREL
  120         FORMAT(7X,'.',A2,'.')
              CALL FMPRNT(MB)
          ELSE
              CALL FMNTRJ(MA,NDIG)
              WRITE (KW,120) LREL
              CALL FMNTRJ(MB,NDIG)
          ENDIF
      ENDIF

!             JCOMP will be 1 if MA.GT.MB
!                           2 if MA.EQ.MB
!                           3 if MA.LT.MB

!             Check for special cases.

      IF (LREL.NE.'EQ' .AND. LREL.NE.'NE' .AND. LREL.NE.'LT' .AND.  &
          LREL.NE.'GT' .AND. LREL.NE.'LE' .AND. LREL.NE.'GE') THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          IF (KWARN.LE.0) GO TO 170
          WRITE (KW,130) LREL
  130     FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE',  &
                 ' FMCOMP'//1X,A,' IS NOT ONE OF THE SIX RECOGNIZED',  &
                 ' COMPARISONS.'//' .FALSE. HAS BEEN RETURNED.'/)
          GO TO 170
      ENDIF

      IF (MA(1).EQ.KUNKNO .OR. MB(1).EQ.KUNKNO) THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          GO TO 170
      ENDIF

      IF (ABS(MA(1)).EQ.KEXPOV .AND. MA(1).EQ.MB(1) .AND.  &
          MA(2).EQ.MB(2)) THEN
          KFLAG = -4
          FMCOMP = .FALSE.
          IF (KWARN.LE.0) GO TO 170
          WRITE (KW,140)
  140     FORMAT(/' ERROR OF TYPE KFLAG = -4 IN FM PACKAGE IN ROUTINE',  &
                 ' FMCOMP'//' TWO NUMBERS IN THE SAME OVERFLOW OR',  &
                 ' UNDERLOW CATEGORY CANNOT BE COMPARED.'//  &
                 ' .FALSE. HAS BEEN RETURNED.'/)
         GO TO 170
      ENDIF

!             Check for zero.

      KFLAG = 0
      IF (MA(2).EQ.0) THEN
          JCOMP = 2
          IF (MB(2).LT.0) JCOMP = 1
          IF (MB(2).GT.0) JCOMP = 3
          GO TO 160
      ENDIF
      IF (MB(2).EQ.0) THEN
          JCOMP = 1
          IF (MA(2).LT.0) JCOMP = 3
          GO TO 160
      ENDIF
!             Check for opposite signs.

      IF (MA(2).GT.0 .AND. MB(2).LT.0) THEN
          JCOMP = 1
          GO TO 160
      ENDIF
      IF (MB(2).GT.0 .AND. MA(2).LT.0) THEN
          JCOMP = 3
          GO TO 160
      ENDIF

!             See which one is larger in absolute value.

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

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

      JCOMP = 2

!             Now match the JCOMP value to the requested comparison.

  160 IF (JCOMP.EQ.1 .AND. MA(2).LT.0) THEN
          JCOMP = 3
      ELSE IF (JCOMP.EQ.3 .AND. MB(2).LT.0) THEN
          JCOMP = 1
      ENDIF

      FMCOMP = .FALSE.
      IF (JCOMP.EQ.1 .AND. (LREL.EQ.'GT' .OR. LREL.EQ.'GE' .OR.  &
                            LREL.EQ.'NE')) FMCOMP = .TRUE.

      IF (JCOMP.EQ.2 .AND. (LREL.EQ.'EQ' .OR. LREL.EQ.'GE' .OR.  &
                            LREL.EQ.'LE')) FMCOMP = .TRUE.

      IF (JCOMP.EQ.3 .AND. (LREL.EQ.'NE' .OR. LREL.EQ.'LT' .OR.  &
                            LREL.EQ.'LE')) FMCOMP = .TRUE.

  170 IF (NTRACE.NE.0) THEN
          IF (NCALL.LE.LVLTRC .AND. ABS(NTRACE).GE.1) THEN
              IF (KFLAG.EQ.0) THEN
                  WRITE (KW,180) NCALL,JBASE,NDIG
  180             FORMAT(' FMCOMP',15X,'CALL LEVEL =',I2,5X,'JBASE =',  &
                         I10,5X,'NDIG =',I6)
              ELSE
                  WRITE (KW,190) NCALL,JBASE,NDIG,KFLAG
  190             FORMAT(' FMCOMP',6X,'CALL LEVEL =',I2,4X,'JBASE =',  &
                         I10,4X,'NDIG =',I6,4X,'KFLAG =',I3)
              ENDIF
              IF (FMCOMP) THEN
                  WRITE (KW,200)
  200