PROGRAM TEST90
 
!                             David M. Smith              9-17-96
 
!    Program using the FM Fortran-90 modules for doing
!    arithmetic using the FM, IM, and ZM derived types.
 
!    This program does the same calculations as FMSAMPLE and ZMSAMPLE.
 
!    The output is saved in file SAMPLE90.LOG.  A comparison file,
!    SAMPLE90.CHK, is provided showing the expected output from 32-bit
!    (IEEE arithmetic) machines.  When run on other computers, all the
!    numerical results should still be the same, but the number of terms
!    needed for some of the results might be slightly different.  The
!    program checks all the results and the last line of the log file
!    should be "All results were ok."
 
!    In a few places, an explicit call is made to an FM or ZM routine.
!    For a call like CALL FM_FORM('F65.60',MAFM,ST1), note that the
!    "FM_" form is used since MAFM is a TYPE (FM) variable and not just
!    an array.  See the discussion in FMZM90.f.
 
 
      USE FMZM
 
      TYPE ( FM ) MAFM,MBFM,MCFM,MDFM
      TYPE ( IM ) MAIM,MBIM,MCIM
      TYPE ( ZM ) MAZM,MBZM,MCZM,MDZM
 
!             Character string used to format multiple-precision output.
 
      CHARACTER *80 ST1
      INTEGER ITER,J,K,KLOG,NERROR
 
!             Note that any program using the FM package MUST call
!             FM_SET before using the package.
 
!             Since the argument to FM_SET is not an FM number,
!             calling FMSET would have the same effect.  The "FM_"
!             version is available so that all calls in a program
!             using the derived types can have the "FM_" form.
 
!             Later in this program complex arithmetic is be used,
!             and ZM_SET is called there to initialize the ZM package.
 
!             Set precision to give at least 60 significant digits
!             and initialize the FMLIB package.
 
      CALL FM_SET(60)
 
      NERROR = 0
 
!             Write output to the standard FM output (unit KW, defined
!             in subroutine FMSET), and also to the file SAMPLE90.LOG.
 
      KLOG = 18
      OPEN (KLOG,FILE='SAMPLE90.LOG')
 
!             1.  Find a root of the equation
!                 f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0.
 
!                 Use Newton's method with initial guess x = 3.12.
!                 This version is not tuned for speed.  See the FMSQRT
!                 routine for possible ways to increase speed.
!                 Horner's rule is used to evaluate the function.
 
!                 MAFM is the previous iterate.
!                 MBFM is the current iterate.
 
!             TO_FM is a function for converting other types of numbers
!             to type FM.  Note that TO_FM(3.12) converts the REAL
!             constant to FM, but it is accurate only to single
!             precision.  TO_FM(3.12D0) agrees with 3.12 to double
!             precision accuracy, and TO_FM('3.12') or
!             TO_FM(312)/TO_FM(100) agrees to full FM accuracy.
 
      MAFM = TO_FM('3.12')
 
!                 Print the first iteration.
 
      WRITE (KW,110)
      WRITE (KLOG,110)
  110 FORMAT(//' Sample 1.  Real root of f(x) = x**5 - 3x**4 + ', &
               'x**3 - 4x**2 + x - 6 = 0.'/// &
               ' Iteration       Newton Approximation')
      CALL FM_FORM('F65.60',MAFM,ST1)
      WRITE (KW,120) 0,ST1(1:65)
      WRITE (KLOG,120) 0,ST1(1:65)
  120 FORMAT(/I10,4X,A)
 
      DO ITER = 1, 10
 
!                 MCFM is f(MAFM).
 
         MCFM = ((((MAFM-3)*MAFM+1)*MAFM-4)*MAFM+1)*MAFM-6
 
!                 MDFM is f'(MAFM).
 
         MDFM = (((5*MAFM-12)*MAFM+3)*MAFM-8)*MAFM+1
 
         MBFM = MAFM - MCFM/MDFM
 
!                 Print each iteration.
 
         CALL FM_FORM('F65.60',MBFM,ST1)
         WRITE (KW,120) ITER,ST1(1:65)
         WRITE (KLOG,120) ITER,ST1(1:65)
 
!                 Stop iterating if MAFM and MBFM agree to over
!                 60 places.
 
         MDFM = ABS(MAFM-MBFM)
         IF (MDFM.LT.1.0D-61) EXIT
 
!                 Set MAFM = MBFM for the next iteration.
 
         MAFM = MBFM
      ENDDO
 
!                 Check the answer.
 
      MCFM = TO_FM('3.120656215326726500470956013523797484654623'// &
                   '9355990660149888284358')
      IF (ABS(MCFM-MBFM).GT.1.0D-61) THEN
          NERROR = NERROR + 1
          WRITE (KW,130)
          WRITE (KLOG,130)
  130     FORMAT(/' Error in sample case number 1.'/)
      ENDIF
 
!             2.  Compute the Riemann Zeta function for s=3.
 
!                 Use Gosper's formula:  Zeta(3) =
!                 (5/4)*Sum[ (-1)**k * (k!)**2 / ((k+1)**2 * (2k+1)!) ]
!                 while k = 0, 1, ....
 
!                 MAFM is the current partial sum.
!                 MBFM is the current term.
!                 MCFM is k!
!                 MDFM is (2k+1)!
 
      MAFM = 1
      MCFM = 1
      MDFM = 1
      DO K = 1, 200
         MCFM = K*MCFM
         J = 2*K*(2*K+1)
         MDFM = J*MDFM
         MBFM = MCFM**2
         J = (K+1)*(K+1)
         MBFM = (MBFM/J)/MDFM
         IF (MOD(K,2).EQ.0) THEN
             MAFM = MAFM + MBFM
         ELSE
             MAFM = MAFM - MBFM
         ENDIF
 
!                 Test for convergence.
 
         IF (MAFM-MBFM.EQ.MAFM) THEN
             WRITE (KW,140) K
             WRITE (KLOG,140) K
  140        FORMAT(///' Sample 2.',8X,I5,' terms were added'/)
             EXIT
         ENDIF
      ENDDO
 
!                 Print the result.
 
      MAFM = (5*MAFM)/4
      CALL FM_FORM('F65.60',MAFM,ST1)
      WRITE (KW,150) ST1(1:65)
      WRITE (KLOG,150) ST1(1:65)
  150 FORMAT(' Zeta(3) = ',A)
 
!                 Check the answer.
 
      MCFM = TO_FM('1.20205690315959428539973816151144999076498'// &
                   '6292340498881792271555')
      IF (ABS(MAFM-MCFM).GT.1.0D-61) THEN
          NERROR = NERROR + 1
          WRITE (KW,160)
          WRITE (KLOG,160)
  160     FORMAT(/' Error in sample case number 2.'/)
      ENDIF
 
!             3.  Integer multiple precision calculations.
 
!                 Fermat's theorem says  x**(p-1) mod p = 1
!                 when p is prime and x is not a multiple of p.
!                 If  x**(p-1) mod p  gives 1 for some p with
!                 several different x's, then it is very likely
!                 that p is prime (but it is not certain until
!                 further tests are done).
 
!                 Find a 70-digit number p that is "probably" prime.
 
!                 MAIM is the value p being tested.
 
      MAIM = TO_IM(10)**69
 
!                 To speed up the search, test only values that are
!                 not multiples of 2, 3, 5, 7, 11, 13.
 
      K = 2*3*5*7*11*13
      MAIM = (MAIM/K)*K + K + 1
      MCIM = 3
 
      DO J = 1, 100
 
!                 Compute 3**(p-1) mod p
 
         MBIM = MAIM - 1
         CALL IM_PMOD(MCIM,MBIM,MAIM,MCIM)
         IF (MCIM.EQ.1) THEN
 
!                 Check that 7**(p-1) mod p is also 1.
 
             MCIM = 7
             CALL IM_PMOD(MCIM,MBIM,MAIM,MCIM)
             IF (MCIM.EQ.1) THEN
                 WRITE (KW,170) J
                 WRITE (KLOG,170) J
  170            FORMAT(///' Sample 3.',8X,I5,' values were tested'/)
                 EXIT
             ENDIF
         ENDIF
 
         MCIM = 3
         MAIM = MAIM + K
      ENDDO
 
!                 Print the result.
 
      CALL IM_FORM('I72',MAIM,ST1)
      WRITE (KW,180) ST1(1:72)
      WRITE (KLOG,180) ST1(1:72)
  180 FORMAT(' p = ',A)
 
!                 Check the answer.
 
      MCIM = TO_IM('1000000000000000000000000000000000000000000'// &
                   '000000000000000000000659661')
      IF (MAIM.NE.MCIM) THEN
          NERROR = NERROR + 1
          WRITE (KW,190)
          WRITE (KLOG,190)
  190     FORMAT(/' Error in sample case number 3.'/)
      ENDIF
 
 
!             Complex arithmetic.
 
!             Set precision to give at least 30 significant digits
!             and initialize the ZMLIB package.  Both FM and ZM
!             operations will now have this precision.
 
!             Note that any program using the ZM package MUST call
!             ZM_SET before using the package.
 
      CALL ZM_SET(30)
 
!             4.  Find a complex root of the equation
!                 f(x) = x**5 - 3x**4 + x**3 - 4x**2 + x - 6 = 0.
 
!                 Newton's method with initial guess x = .56 + 1.06 i.
!                 This version is not tuned for speed.  See the ZMSQRT
!                 routine for possible ways to increase speed.
!                 Horner's rule is used to evaluate the function.
 
!                 MAZM is the previous iterate.
!                 MBZM is the current iterate.
 
      MAZM = TO_ZM('.56 + 1.06 i')
 
!                 Print the first iteration.
 
      WRITE (KW,200)
      WRITE (KLOG,200)
  200 FORMAT(//' Sample 4.  Complex root of f(x) = x**5 - 3x**4 + ', &
               'x**3 - 4x**2 + x - 6 = 0.'/// &
               ' Iteration       Newton Approximation')
      CALL ZM_FORM('F32.30','F32.30',MAZM,ST1)
      WRITE (KW,210) 0,ST1(1:69)
      WRITE (KLOG,210) 0,ST1(1:69)
  210 FORMAT(/I6,4X,A)
 
      DO ITER = 1, 10
 
!                 MCZM is f(MAZM).
 
         MCZM = ((((MAZM-3)*MAZM+1)*MAZM-4)*MAZM+1)*MAZM-6
 
!                 MDZM is f'(MAZM).
 
         MDZM = (((5*MAZM-12)*MAZM+3)*MAZM-8)*MAZM+1
 
         MBZM = MAZM - MCZM/MDZM
 
!                 Print each iteration.
 
         CALL ZM_FORM('F32.30','F32.30',MBZM,ST1)
         WRITE (KW,210) ITER,ST1(1:69)
         WRITE (KLOG,210) ITER,ST1(1:69)
 
!                 Stop iterating if MAZM and MBZM agree to over
!                 30 places.
 
         IF (ABS(MAZM-MBZM).LT.1.0D-31) EXIT
 
!                 Set MAZM = MBZM for the next iteration.
 
         MAZM = MBZM
      ENDDO
 
!                 Check the answer.
 
      MCZM = TO_ZM('0.561958308335403235498111195347453 +'// &
                   '1.061134679604332556983391239058885 i')
      IF (ABS(MCZM-MBZM).GT.1.0D-31) THEN
          NERROR = NERROR + 1
          WRITE (KW,220)
          WRITE (KLOG,220)
  220     FORMAT(/' Error in sample case number 4.'/)
      ENDIF
 
!             5.  Compute Exp(1.23-2.34i).
 
!                 Use the direct Taylor series.  See the ZMEXP routine
!                 for a faster way to get Exp(x).
 
!                 MAZM is x.
!                 MBZM is the current term, x**n/n!.
!                 MCZM is the current partial sum.
 
      MAZM = TO_ZM('1.23-2.34i')
      MBZM = 1
      MCZM = 1
      DO K = 1, 100
         MBZM = MBZM*MAZM/K
         MDZM = MCZM + MBZM
 
!                 Test for convergence.
 
         IF (MDZM.EQ.MCZM) THEN
             WRITE (KW,230) K
             WRITE (KLOG,230) K
  230        FORMAT(///' Sample 5.',8X,I5,' terms were added ', &
                    'to get Exp(1.23-2.34i)'/)
             EXIT
         ENDIF
         MCZM = MDZM
      ENDDO
 
!                 Print the result.
 
      CALL ZM_FORM('F33.30','F32.30',MCZM,ST1)
      WRITE (KW,240) ST1(1:70)
      WRITE (KLOG,240) ST1(1:70)
  240 FORMAT(' Result= ',A)
 
!                 Check the answer.
 
      MDZM = TO_ZM('-2.379681796854777515745457977696745 -'// &
                   ' 2.458032970832342652397461908326042 i')
      IF (ABS(MDZM-MCZM).GT.1.0D-31) THEN
          NERROR = NERROR + 1
          WRITE (KW,250)
          WRITE (KLOG,250)
  250     FORMAT(/' Error in sample case number 5.'/)
      ENDIF
 
      IF (NERROR.EQ.0) THEN
          WRITE (KW,260) ' All results were ok.'
          WRITE (KLOG,260) ' All results were ok.'
  260     FORMAT(//A/)
      ENDIF
      END